///////////////////////////////////////////////////////////////////////////
//																																			 //
//	Class: Enriched_polyhedron
//																																			 //
///////////////////////////////////////////////////////////////////////////

#pragma once
#pragma  warning(disable : 4273)
#pragma  warning(disable : 4503)
#pragma  warning(disable : 4018)
// CGAL	stuff
// #pragma comment(lib,"glew32s.lib")
// //#pragma comment(lib,"glew32s.lib")
#include <CGAL/Cartesian.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/HalfedgeDS_vertex_max_base_with_id.h>
#include <CGAL/HalfedgeDS_halfedge_max_base_with_id.h>
#include <CGAL/HalfedgeDS_face_max_base_with_id.h>
#include <CGAL/HalfedgeDS_decorator.h> 
#include <CGAL/Random_access_adaptor.h>
#include <CGAL/intersections.h>
#include <CGAL/Kd_tree.h>
#include <CGAL/Fuzzy_sphere.h>
#include <CGAL/Search_traits_3.h>
#include <list>
#include <set>
#include <stack>
#include <map>
#include <queue>
#include <string>
#include "Math\mathlib.h"
#include "Math\matrix.h"
#include "TriangleMesh\curvature.h"
#include "TriangleMesh\file_loader\builder.h"
#include "GLAddon\Utility.h"
#include "GLAddon\PsudoColorRGB.h"
#include <fstream>
#include "GLAddon\OpenGLFunc.h"

// linear solver
//#include <CGAL/Cholmod_matrix.h>
//#include <CGAL/Cholmod_vector.h>

#define VH(v) index_to_vertex_map[(v)]					/*get V handle*/
#define VP(v) index_to_vertex_map[(v)]->point()			/*get V point*/
#define VN(v) index_to_vertex_map[(v)]->normal()		/*get V normal*/
#define VB(v) index_to_vertex_map[(v)]->vertex_begin()  /*get the 1st V handle*/
#define FH(f) index_to_facet_map[(f)]					/*get F handle*/
#define FN(f) index_to_facet_map[(f)]->normal()         /*get F normal*/
#define FB(f) index_to_facet_map[(f)]->facet_begin()    /*get the 1st F handle*/
#define EH(e) index_to_halfedge_map[(e)]			    /*get halfedge handle*/

/*convert cgal point to vector3*/
#define P2V3(p)	  Vector3( (p).x(), (p).y(), (p).z()) 
#define V32P(v)	  Point( (v).x, (v).y, (v).z) 
//	return vct;}\						

#define BN (1.0e8)
#define SQBN (1.0e4)
#define PICKED 4      //for add info bit to the tag
// a refined facet with a normal and a tag
template <class	Refs, class T, class P,	class Norm>
class Enriched_facet : public CGAL::HalfedgeDS_face_max_base_with_id<Refs, CGAL::Tag_false, std::size_t>
{
	// tag
	int	m_tag;
	// normal
	Norm m_normal;

public:
	// life	cycle
	// no constructors to repeat, since only default constructor mandatory
	Enriched_facet() {}

	// tag
	const int& tag() const {return m_tag;}
	int& tag() {return m_tag;}
	void tag(const int& t)	{m_tag = t;}
	void add_tag(const unsigned int& t) { m_tag |= (1<<t); }
	void remove_tag(const unsigned int& t) { m_tag &= ~(1<<t); }
	int have_tag(const unsigned int& t) { return (m_tag & (1<<t)); }
	// normal
	typedef	Norm Normal_3;
	Normal_3& normal() { return	m_normal; }
	const Normal_3& normal() const { return m_normal; }
};

// a refined halfedge with a general tag and 
// a binary tag to indicate wether it belongs 
// to the control mesh or not
template <class	Refs, class Tprev, class Tvertex, class Tface, class Norm>
class Enriched_halfedge : public CGAL::HalfedgeDS_halfedge_max_base_with_id<Refs, std::size_t>
{
private:

	// tag
	int	m_tag; 

	// option	for	edge superimposing
	bool m_control_edge;

	// normal
	Norm m_vnormal;

public:

	// life	cycle
	Enriched_halfedge()
	{
		m_control_edge = true;
	}

	// tag
	const int& tag() const {return m_tag;}
	int& tag() {return m_tag;}
	void tag(const int& t)	{m_tag = t;}
	void add_tag(const unsigned int& t) { m_tag |= (1<<t); }
	void remove_tag(const unsigned int& t) { m_tag &= ~(1<<t); }
	int have_tag(const unsigned int& t) { return (m_tag & (1<<t)); }
	// control edge	
	bool& control_edge() {return m_control_edge;}
	const bool& control_edge() const {return m_control_edge;}
	void control_edge(const bool flag) {m_control_edge = flag;}

	// normal
	Norm& vnormal() { return m_vnormal; }
	const Norm& vnormal() const { return m_vnormal; }
};

// a refined vertex with a normal and a tag
template <class	Refs, class T, class P, class N, class Kernel>
class Enriched_vertex : public CGAL::HalfedgeDS_vertex_max_base_with_id<Refs, P, std::size_t>
{
public:
	typedef	typename Kernel::FT FT;
	typedef	typename N Normal;
	typedef typename CCurvature<Kernel> Curvature;

private:
	int	m_tag;
	short m_dup;	// need duplicating? If a vertex has more than one attributes, it must
	//	be duplicated in building vertex buffer for VBO rendering
	Normal m_normal;
	Curvature m_curvature;
	union 
	{
		struct { FT m_u, m_v; };
		struct { FT m_texcoord[2]; };
	};	// texture coordinates u, v

public:
	// life	cycle
	Enriched_vertex() {}
	// repeat	mandatory	constructors
	Enriched_vertex(const P& pt)
		:	CGAL::HalfedgeDS_vertex_max_base_with_id<Refs, P, std::size_t>(pt)
	{
	}

	// normal
	Normal& normal() { return m_normal; }
	const Normal& normal() const { return m_normal; }

	// curvature
	Curvature& curvature() { return m_curvature; }
	const Curvature& curvature() const { return m_curvature; }

	// tag
	int& tag() { return m_tag; }
	const int& tag() const { return m_tag; }
	void tag(const int& t)	{		m_tag	= t; }
	void add_tag(const unsigned int& t) { 
		m_tag |= (1<<t); }
	void remove_tag(const unsigned int& t) { 
		m_tag &= ~(1<<t); }
	int have_tag(const unsigned int& t) { return (m_tag & (1<<t)); }
	// dup
	short& dup() { return m_dup; }
	const short& dup() const { return m_dup; }
	void dup(const short& d)	{ m_dup	= d; }

	// texcoord
	FT& u() { return m_u; }
	FT& v() { return m_v; }
	FT* texcoord() { return &m_texcoord[0]; }
	const FT& u() const { return m_u; }
	const FT& v() const { return m_v; }
	const FT* texcoord() const { return &m_texcoord[0]; }
	void u(const FT uc) { m_u = uc; }
	void v(const FT vc) { m_v = vc; }
	void texcoord(const FT uc, const FT vc) { m_u = uc; m_v = vc; }
	void texcoord(const FT* const tc) { m_u = tc[0]; m_v = tc[1]; }
};

// A redefined items class for the Polyhedron_3	
// with	a refined vertex class that	contains a 
// member for the normal vector	and	a refined
// facet with a normal vector instead of the 
// plane equation (this is an alternative	
// solution	instead	of using 
// Polyhedron_traits_with_normals_3).

struct Enriched_items : public CGAL::Polyhedron_items_3
{
	// wrap	vertex
	template <class Refs, class Traits>
	struct Vertex_wrapper
	{
		typedef	typename Traits::Point_3	Point;
		typedef	typename Traits::Vector_3	Normal;
		typedef	Enriched_vertex<Refs, CGAL::Tag_true, Point, Normal, Traits> Vertex;
	};

	// wrap	face
	template <class Refs, class Traits>
	struct Face_wrapper
	{
		typedef	typename Traits::Point_3	Point;
		typedef	typename Traits::Vector_3	Normal;
		typedef	Enriched_facet<Refs, CGAL::Tag_true, Point, Normal> Face;
	};

	// wrap	halfedge
	template <class Refs, class Traits>
	struct Halfedge_wrapper
	{
		typedef	typename Traits::Vector_3	Normal;
		typedef	Enriched_halfedge<Refs, CGAL::Tag_true, CGAL::Tag_true, CGAL::Tag_true, Normal>	Halfedge;
	};
};

template <class	kernel, class items>
class ROI_vertex;

//*********************************************************
template <class	kernel, 
class items, 
#ifndef CGAL_CFG_NO_TMPL_IN_TMPL_PARAM
	template <class T, class I, class A>
#endif
class T_HDS = CGAL::HalfedgeDS_default,
class A = CGAL_ALLOCATOR(int)>
class Enriched_polyhedron : public CGAL::Polyhedron_3<kernel,items>
{
public :
	typedef	typename kernel::FT			FT;
	typedef	typename kernel::Point_3	Point;
	typedef	typename kernel::Vector_3	Vector;
	typedef typename Vertex::Normal		Normal;
	typedef CGAL::Random_access_adaptor<Vertex_iterator>	Random_access_vertex_index;
	typedef CGAL::Random_access_adaptor<Facet_iterator>		Random_access_facet_index;
	typedef CGAL::Random_access_adaptor<Halfedge_iterator>	Random_access_halfedge_index;
	typedef CGAL::HalfedgeDS_decorator<HalfedgeDS>			HalfedgeDS_decorator;
	typedef std::stack<std::pair<int, std::list<int>>>		Stack_list;
	typedef CGAL::Search_traits_3<kernel>					STraits;
	typedef CGAL::Kd_tree<STraits>							KD_tree;
	typedef CGAL::Fuzzy_sphere<STraits>						Fuzzy_sphere;
#ifdef CGAL_USE_TAUCS
	//	typedef CGAL::Taucs_vector<FT>				TVector;
	//	typedef CGAL::Taucs_matrix<FT>				TMatrix;
	//	typedef CGAL::Taucs_symmetric_matrix<FT>	TSMatrix;
#endif
//	typedef CGAL::Cholmod_matrix<double>			CMatrix;	// CHOLMOD only supports double floating-point number
//	typedef CGAL::Cholmod_vector<double>			CVector;	// CHOLMOD only supports double floating-point number
	typedef std::vector<std::set<int>>				VEC_SI;
	typedef std::vector<std::list<Facet_handle>>	VEC_LFh;
	typedef std::map<int,std::pair<short,short>>	MAP_ISP;	// map int to short pair
	typedef std::map<int,std::pair<Vector3,short>>	MAP_IVS;	// map int to Vector3&short

public:
	std::list<int> picked_vertices;
	std::list<int> picked_facets;
	std::list<int> picked_edges;	
	CString	m_fileName; 
private:
	// bounding box
	FT m_min[3];
	FT m_max[3];
	FT m_scale;
	FT m_center[3];
	GLuint m_gl_bb_disp_list1;
	GLuint m_gl_bb_disp_list2;
	FT m_suitable_point_size;  //suitalbe size for points on surface

	bool m_is_set_items_id;
	bool m_is_set_id_to_items_map;

	// type
	bool m_pure_quad;
	bool m_pure_triangle;

	GLuint m_gl_display_list;
	GLuint m_gl_vbo_pos;
	GLuint m_gl_vbo_norm;
	GLuint m_gl_vbo_texcoord;
	GLuint m_gl_vbo_color;		// for drawing handles
	GLuint m_gl_vbo_hf_color;
	GLuint m_gl_vbo_hv;			// harmonic values buffer object
	GLuint m_gl_ibo;
	std::vector<GLuint> m_gl_ibo_mf_hdl;
	std::vector<int> m_vb_id;	// vertex id of vertices of vertex buffer
	GLuint *m_ibo_array;
	std::vector<int> m_vb_ia_id;// id (in ibo array) of vertices of vertex buffer

	std::list<Facet_handle>		m_roi_f;			// facets in roi

	FT m_min_mean_curvature;
	FT m_max_mean_curvature;

	// Deformation
	ROI_vertex<kernel,items>	*m_roi_v;
	std::map<int, int>			m_gi_to_ri_map;		// map global index to roi index 
	int							m_nb_roi_v;
	int							m_nb_broi_v;
	int							m_roi_hdl_v;		// roi index of handle
	int							m_roi_hdl_gi;
	FT							m_roi_hdl_r;		// radia
	//
	std::set<int>				m_fixed_hdl_v;		// fixed handle vertices
	std::set<int>				m_manipl_hdl_v;		// manipulating handle vertices
	//////////////////////////////////////////////////////////////////////////
	short						*m_hdl;				// handle config: handle no. of each vertex
	short						*m_hdl_tag;			// handle tag:  0: up-to-date -1: need to be downdated
	//				1: newly added or need to be updated
	//				2: non-constraint
	VEC_SI						m_hdl_v;			// vertex ids in each handle
	VEC_LFh						m_hdl_f;			// facets handle in each handle
	int							m_curr_manip;		// current manipulated handle
	bool						m_is_comp_mh_hf;
	bool						m_is_mh_hf_utd;
	std::vector<FT*>			m_mh_hf;
	//////////////////////////////////////////////////////////////////////////
	std::vector<std::set<int>>	m_joint_v;			// vertices in joint set
	std::list<Facet_handle>		m_fixed_hdl_f;		// fixed handle facets
	std::list<Facet_handle>		m_manipl_hdl_f;		// manipulating handle facets
	int							m_pick_hdl_f;		// picked handle facet
	Vector3						m_pick_hdl_pos;
	Vector3						*m_dif_coords;
	int							m_d_LM;				// Laplacian matrix dimension
	FT							m_pc_weight;		// weight of positional constraint
	bool						m_is_set_roi;
	bool						m_is_set_roi_hdl;	// have set handle vertex?
	bool						m_is_init_roi_df;
	bool						m_is_comp_roi_hf;	// have computed harmonic field
	bool						m_is_init_hb_df;
	bool						m_is_set_manipl_hdl;
	bool						m_is_set_fixed_hdl;
	bool						m_is_comp_hb_hf;
#ifdef CGAL_USE_TAUCS
	taucs_ccs_matrix			*m_LM;				// Laplacian matrix
	TMatrix						*m_tmAM;			// associated matrix
	TVector						*m_tvROIDC[3];		// roi differential coordinates
	TVector						*m_tvB;				// B vector
	TVector						*m_tvX[3];			// X vector
#endif
	void						*m_CF;				// Cholesky factorization
	//cholmod_factor				*m_L;				// factor
	FT							*m_roi_hf;			// harmonic field in roi
	FT							*m_hb_hf;			// handle-based hf over the entire mesh
	bool						m_is_unitized;
	FT							m_unitize_sf;
	MAP_ISP						m_hf_sites;			// flag: 1: source  0: sink
	// tag: -1: downdate  0: up-to-date  1: update
	int							m_pick_hf_site;
	FT							*m_s_hf;			// site-based harmonic field
	//cholmod_factor				*m_s_L;				// factor
	bool						m_is_built_hf_s_vbo;
	bool						m_is_utd_sb_hf;
	float						m_hf_site_radii;

	MAP_IVS						m_hvf_sites;		// harmonic vector fields, tag: -1:downdate 0:up-to-date 1:update
	int							m_pick_hvf_site;
	bool						m_is_built_hvf_s_vbo;
	bool						m_is_utd_sb_hvf;
	bool						m_is_comp_sb_hvf;
	Vector3						*m_sb_vf;			// site-based harmonic vector field
	FT							m_s_hvm_max;		// maximum magnitude
	FT							m_s_hvm_min;		// minimum magnitude
	//cholmod_factor				*m_s_hvf_L;			// factor
	float						m_hvf_site_radii;
	std::vector<int>			m_s_hvf_p;			// sampled points for illustrating site-based harmonic vector field

	std::list<int>				m_curve;
	std::list<int>				m_c_add, m_c_rm;
	int							m_cd_head, m_cd_curr, m_cd_tail; // draw head and draw tail
	short						m_c_status;
	int							m_seed;
	std::vector<int>			m_hole_v;
	//cholmod_factor				*m_mbr_F;			// factor for computing membrane
	//cholmod_factor				*m_mbr_F_e;			// factor for evolving membrane
	//cholmod_sparse				*m_mbr_Lt;			// L^t for evolving membrane
	bool						m_is_mbr_utd;
	Vector3						*m_mbr_dc;

	FT							m_avg_edge_len;

	GLubyte						*m_texture_mat;		// material per vertex for synthesized texture
	GLuint						m_gl_vbo_tex_color;
	bool						m_is_built_tex_vbo;
	std::vector<Vector3>		m_points_copy;
	bool						m_bPointsCopied;
public:
	Random_access_vertex_index	index_to_vertex_map;
	Random_access_facet_index	index_to_facet_map;
	Random_access_halfedge_index    index_to_halfedge_map;
	// life	cycle
	Enriched_polyhedron()	
	{
		m_pure_quad = false;
		m_pure_triangle = false;
		m_roi_v = NULL;
		picked_vertices.clear();
		picked_facets.clear();
		picked_edges.clear();
		m_dif_coords = NULL;
		m_d_LM = NULL;
		m_pc_weight = 100.0f;
		m_is_set_roi = false;
		m_is_set_roi_hdl = false;
		m_fixed_hdl_v.clear();
		m_manipl_hdl_v.clear();
		m_hdl_v.clear();
		m_hdl_f.clear();
		m_hdl = NULL;
		m_hdl_tag = NULL;
		m_curr_manip = -1;
		m_mh_hf.clear();
		m_is_comp_mh_hf = false;
		m_is_mh_hf_utd = false;
		//m_L = NULL;
		m_roi_hdl_v = -1;
		m_is_init_hb_df = false;
		m_pick_hdl_f = -1;
#ifdef CGAL_USE_TAUCS
		m_LM = NULL;
		m_tmAM = NULL;
		m_tvROIDC[0] = m_tvROIDC[1] = m_tvROIDC[2] = NULL;
		m_tvB = NULL;
		m_tvX[0] = m_tvX[1] = m_tvX[2] = NULL;
		taucs_logfile("taucslog.txt");
#endif
		m_CF = NULL;
		m_roi_hf = NULL;
		m_hb_hf = NULL;
		m_roi_hdl_gi = -1;
		m_roi_hdl_r = 0.25;
		m_is_comp_roi_hf = false;
		m_is_init_roi_df = false;
		m_is_set_items_id = false;
		m_is_set_manipl_hdl = false;
		m_is_set_fixed_hdl = false;
		m_is_comp_hb_hf = false;
		m_is_set_id_to_items_map = false;
		m_is_unitized = false;
		m_unitize_sf = 5.0;
		m_gl_vbo_pos = 0;
		m_gl_vbo_norm = 0;
		m_gl_vbo_hf_color = 0;
		m_gl_vbo_tex_color = 0;
		m_gl_vbo_color = 0;
		m_gl_vbo_hv = 0;
		m_gl_ibo = 0;
		m_ibo_array = NULL;
		m_hf_sites.clear();
		m_s_hf = NULL;
		//m_s_L = NULL;
		m_is_built_hf_s_vbo = false;
		m_is_built_tex_vbo = false;
		m_is_utd_sb_hf = false;
		m_hvf_sites.clear();
		m_pick_hf_site = 0;
		m_pick_hvf_site = 0;
		m_is_built_hvf_s_vbo = false;
		m_is_utd_sb_hvf = true;
		m_is_comp_sb_hvf = false;
		m_s_hvf_p.clear();
		m_sb_vf = NULL;
		//m_s_hvf_L = NULL;
		m_hf_site_radii = 0.06;
		m_hvf_site_radii = 0.03;
		m_curve.clear(); m_c_add.clear(); m_c_rm.clear();
		m_cd_head = m_cd_curr = m_cd_tail = 0;
		m_c_status = 0;
		m_seed = -1;
		m_hole_v.clear();
		//m_mbr_F = NULL; m_mbr_F_e = NULL;
		m_is_mbr_utd = false;
		m_mbr_dc = NULL;
		//m_mbr_Lt = NULL;
		m_avg_edge_len = 0.0;
		m_texture_mat = NULL;
		m_s_hvm_min = 0.0;
		m_s_hvm_max = 0.0;
		m_bPointsCopied = false;
	}
	Enriched_polyhedron(const Enriched_polyhedron& poly)
		: Polyhedron_3(poly)
	{
		m_pure_quad = false;
		m_pure_triangle = false;
		m_roi_v = NULL;
		picked_vertices.clear();
		picked_facets.clear();
		picked_edges.clear();
		m_dif_coords = NULL;
		m_d_LM = NULL;
		m_pc_weight = 100.0f;
		m_is_set_roi_hdl = false;
		m_is_set_roi = false;
		m_L = NULL;
		m_fixed_hdl_v.clear();
		m_manipl_hdl_v.clear();
		m_hdl = NULL;
		m_hdl_tag = NULL;
		m_hdl_v.clear();
		m_hdl_f.clear();
		m_mh_hf.clear();
		m_curr_manip = -1;
		m_is_comp_mh_hf = false;
		m_is_mh_hf_utd = false;
		m_roi_hdl_v = -1;
		m_is_init_hb_df = false;
		m_pick_hdl_f = -1;
#ifdef CGAL_USE_TAUCS
		m_LM = NULL;
		m_tmAM = NULL;
		m_tvROIDC[0] = m_tvROIDC[1] = m_tvROIDC[2] = NULL;
		m_tvB = NULL;
		m_tvX[0] = m_tvX[1] = m_tvX[2] = NULL;
		//	taucs_logfile("taucslog.txt");
#endif
		m_CF = NULL;
		m_roi_hf = NULL;
		m_hb_hf = NULL;
		m_roi_hdl_gi = -1;
		m_roi_hdl_r = 0.25;
		m_is_comp_roi_hf = false;
		m_is_init_roi_df = false;
		m_is_set_items_id = false;
		m_is_set_manipl_hdl = false;
		m_is_set_fixed_hdl = false;
		m_is_comp_hb_hf = false;
		m_is_set_id_to_items_map = false;
		m_is_unitized = false;
		m_unitize_sf = 5.0;
		m_gl_vbo_pos = 0;
		m_gl_vbo_norm = 0;
		m_gl_vbo_hf_color = 0;
		m_gl_vbo_tex_color = 0;
		m_gl_vbo_color = 0;
		m_gl_vbo_hv = 0;
		m_gl_ibo = 0;
		m_ibo_array = NULL;
		m_hf_sites.clear();
		m_s_hf = NULL;
		m_s_L = NULL;
		m_is_built_hf_s_vbo = false;
		m_is_built_tex_vbo = false;
		m_is_utd_sb_hf = false;
		m_hvf_sites.clear();
		m_pick_hf_site = 0;
		m_pick_hvf_site = 0;
		m_is_built_hvf_s_vbo = false;
		m_is_utd_sb_hvf = true;
		m_is_comp_sb_hvf = false;
		m_s_hvf_p.clear();
		m_sb_vf = NULL;
		m_s_hvf_L = NULL;
		m_hf_site_radii = 0.06;
		m_hvf_site_radii = 0.02;
		m_curve.clear(); m_c_add.clear(); m_c_rm.clear();
		m_cd_head = m_cd_curr = m_cd_tail = 0;
		m_c_status = 0;
		m_seed = -1;
		m_hole_v.clear();
		m_mbr_F = NULL; m_mbr_F_e = NULL;
		m_is_mbr_utd = false;
		m_mbr_dc = NULL; m_mbr_Lt = NULL;
		m_avg_edge_len = 0.0;
		m_texture_mat = NULL;
		m_s_hvm_min = 0.0;
		m_s_hvm_max = 0.0;
		vbo_line_draw_setting(1.3, 78,78,78);
	}
	virtual	~Enriched_polyhedron() 
	{
		if (::glIsList(m_gl_display_list)) {
			::glDeleteLists(m_gl_display_list, 1);
		}
		if (m_gl_vbo_pos != 0) {
			glDeleteBuffers(1, &m_gl_vbo_pos);
		}
		if (m_gl_vbo_norm != 0) {
			glDeleteBuffers(1, &m_gl_vbo_norm);
		}
		if (m_gl_vbo_hf_color != 0) {
			glDeleteBuffers(1, &m_gl_vbo_hf_color);
		}
		if (m_gl_vbo_tex_color != 0) {
			glDeleteBuffers(1, &m_gl_vbo_tex_color);
		}
		if (m_gl_vbo_color != 0) {
			glDeleteBuffers(1, &m_gl_vbo_color);
		}
		if (m_gl_vbo_hv != 0) {
			glDeleteBuffers(1, &m_gl_vbo_hv);
		}
		if (m_gl_ibo != 0) {
			glDeleteBuffers(1, &m_gl_ibo);
		}
		for (int i=0; i<static_cast<int>(m_mh_hf.size()) ; i++) {
			delete [] m_mh_hf[i];
		}
#ifdef CGAL_USE_TAUCS
		delete [] m_CF;
		delete m_tvROIDC[0];
		delete m_tvROIDC[1];
		delete m_tvROIDC[2];
		delete m_tvB;
		delete m_tvX[0];
		delete m_tvX[1];
		delete m_tvX[2];
		delete m_tmAM;
		taucs_ccs_free(m_LM);
#endif
// 		cholmod_common c;
// 		cholmod_start(&c);	// start CHOLMOD
// 		cholmod_free_factor(&m_L, &c);
// 		cholmod_free_factor(&m_s_L, &c);
// 		cholmod_free_factor(&m_s_hvf_L, &c);
// 		cholmod_free_factor(&m_mbr_F, &c);
// 		cholmod_free_factor(&m_mbr_F_e, &c);
// 		cholmod_free_sparse(&m_mbr_Lt, &c);
// 		cholmod_finish(&c);
		SAFE_DELETE_ARRAY(m_texture_mat);
		SAFE_DELETE_ARRAY(m_mbr_dc);
		SAFE_DELETE_ARRAY(m_sb_vf);
		SAFE_DELETE_ARRAY(m_s_hf);
		SAFE_DELETE_ARRAY(m_roi_hf);
		SAFE_DELETE_ARRAY(m_hb_hf);
		SAFE_DELETE_ARRAY(m_dif_coords);
		SAFE_DELETE_ARRAY(m_roi_v);
		SAFE_DELETE_ARRAY(m_ibo_array);
		SAFE_DELETE_ARRAY(m_hdl);
	}

	void set_tex_per_vert_mat(float **color)
	{
		if (color == NULL) {
			return;
		}
		if (m_texture_mat == NULL) {
			m_texture_mat = new GLubyte[3*size_of_vertices()];
		}
		for (int i=0; i<size_of_vertices(); i++)
		{
			m_texture_mat[3*i+0] = (GLubyte)color[i][0];
			m_texture_mat[3*i+1] = (GLubyte)color[i][1];
			m_texture_mat[3*i+2] = (GLubyte)color[i][2];
		}
		m_is_built_tex_vbo = false;
	}

	bool is_built_tex_vbo(void)
	{
		return m_is_built_tex_vbo;
	}
	FT mean_curvature_min() { return m_min_mean_curvature; }
	void mean_curvature_min(FT min) { m_min_mean_curvature = min; }
	FT mean_curvature_max() { return m_max_mean_curvature; }
	void mean_curvature_max(FT max) { m_max_mean_curvature = max; }

	bool is_set_items_id() {return m_is_set_items_id;}
	bool is_set_id_to_items_map() {return m_is_set_id_to_items_map;}

	// type
	bool is_pure_triangle() { return m_pure_triangle; }
	bool is_pure_quad() { return m_pure_quad; }

	// normals (per	facet, then	per	vertex)
	void compute_normals_per_facet()
	{
		std::for_each(facets_begin(), facets_end(), Facet_normal2<kernel, Facet>());
	}
	void compute_normals_per_vertex()
	{
		std::for_each(vertices_begin(), vertices_end(), Vertex_normal<kernel, Vertex>());
	}
	void compute_normals_per_he_vertex()
	{
		std::for_each(vertices_begin(),
			vertices_end(),
			HalfedgeVertex_normal<kernel, Vertex>(85.0));
	}
	void compute_normals()
	{
		compute_normals_per_facet();
		compute_normals_per_vertex();
		compute_normals_per_he_vertex();
	}
	std::vector<FT> m_edge_Len;
	void compute_all_edge_len()
	{
		m_edge_Len.clear(); m_edge_Len.resize(size_of_halfedges()/2);
		typedef CGAL::N_step_adaptor< Halfedge_iterator, 2> Edge_iterator;
		int eid = 0;
		for ( Edge_iterator e = halfedges_begin(); e != halfedges_end(); ++e, ++eid){
			const Point& p1 = e->vertex()->point();
			const Point& p2 = e->opposite()->vertex()->point();
			Vector edge = p1 - p2;
			FT edge_len = (FT)std::sqrt(edge*edge);
			m_edge_Len[eid] = edge_len;
		}
	}
	FT get_avg_edge_len(bool bRecalc=false)
	{
		if (m_avg_edge_len!=0.0 && !bRecalc) {
			return m_avg_edge_len;
		}
		Vector edge;
		FT len_edge = 0.0;
		int num_edge = 0;
		Halfedge_iterator pHalfEdge = halfedges_begin();
		CGAL_For_all(pHalfEdge, halfedges_end())
		{
			const Point& p1 = pHalfEdge->vertex()->point();
			const Point& p2 = pHalfEdge->opposite()->vertex()->point();
			edge = p1 - p2;
			len_edge += (FT)std::sqrt(edge*edge);
			num_edge++;
		}
		m_avg_edge_len = len_edge/(FT)num_edge;
		return m_avg_edge_len;
	}
	void reverse_winding()
	{
		inside_out();
		compute_normals();
	}

	void reverse_picked_facets()
	{
	}
	// compute bounding	box
	void compute_bounding_box(void)
	{
		if (size_of_vertices() == 0) {
			return;
		}
		Vertex_iterator pVertex=vertices_begin();
		m_min[0] = pVertex->point().x();
		m_max[0] = pVertex->point().x();
		m_min[1] = pVertex->point().y();
		m_max[1] = pVertex->point().y();
		m_min[2] = pVertex->point().z();
		m_max[2] = pVertex->point().z();
		for (; pVertex!=vertices_end(); pVertex++)
		{
			const Point& p = pVertex->point();
			m_min[0] = std::min(m_min[0], p.x());
			m_min[1] = std::min(m_min[1], p.y());
			m_min[2] = std::min(m_min[2], p.z());
			m_max[0] = std::max(m_max[0], p.x());
			m_max[1] = std::max(m_max[1], p.y());
			m_max[2] = std::max(m_max[2], p.z());
		}
	}

	// bounding box
	FT xmin() { return m_min[0]; }
	FT xmax() { return m_max[0]; }
	FT ymin() { return m_min[1]; }
	FT ymax() { return m_max[1]; }
	FT zmin() { return m_min[2]; }
	FT zmax() { return m_max[2]; }

	void unitize(FT given_sf=0.0)
	{
		compute_bounding_box();
		// calculate center of the model
		m_center[0] = (m_max[0] + m_min[0])*0.5;
		m_center[1] = (m_max[1] + m_min[1])*0.5;
		m_center[2] = (m_max[2] + m_min[2])*0.5;
		if(given_sf==0.0){  //to be scaled by mesh size
			FT w, h, d;
			// calculate model width, height, and depth
			w = m_max[0] - m_min[0];
			h = m_max[1] - m_min[1];
			d = m_max[2] - m_min[2];
			if(w<0||h<0||d<0)
				MSG_BOX_ERROR("BOUNDING BOX SIZE ERROR");
			// calculate unitizing scale factor
			m_scale = m_unitize_sf / std::max(std::max(w, h), d);
		}
		else
			m_scale = given_sf;

		// translate around center then scale
		for (Vertex_iterator pVertex=vertices_begin(); pVertex!=vertices_end(); pVertex++)
		{
			Point& point = pVertex->point();
			point = Point((point.x()-m_center[0])*m_scale,
				(point.y()-m_center[1])*m_scale,
				(point.z()-m_center[2])*m_scale);
		}
		for (int i=0; i<3; i++)
		{
			m_min[i] = (m_min[i] - m_center[i])* m_scale;
			m_max[i] = (m_max[i] - m_center[i])* m_scale;
		}
		m_is_unitized = true;
	}
	void ununitize(void)
	{
		FT scl = 1.0 / m_scale;
		for (Vertex_iterator pVertex=vertices_begin(); pVertex!=vertices_end(); pVertex++)
		{
			Point& point = pVertex->point();
			point = Point(point.x()*scl+m_center[0],
				point.y()*scl+m_center[1],
				point.z()*scl+m_center[2]);
		}
	}
	void apply_rotation(Matrix3 rot)
	{
		Vector3 pos;
		for (Vertex_iterator pVertex=vertices_begin(); pVertex!=vertices_end(); pVertex++)
		{
			Point& point = pVertex->point();
			pos = rot.mul(Vector3(point.x(), point.y(), point.z()));
			point = Point(pos[0], pos[1], pos[2]);
		}
	}
	bool is_unitized(void)
	{
		return m_is_unitized;
	}
	FT	get_unitizing_scale_factor(void)
	{
		return m_unitize_sf; //mesh bounding box's size, = x_max * 2.
	}
	FT	get_scale_factor(void)
	{
		return m_scale;
	}
	const FT* get_origin_center(void)
	{
		return m_center;
	}
	void gl_build_bounding_box_disp_list1(void) //DRAW BBOX in stipple mode
	{
		if (size_of_vertices() == 0) {
			return;
		}
		FT w, h, d;
		// calculate model width, height, and depth
		w = m_max[0] - m_min[0];
		h = m_max[1] - m_min[1];
		d = m_max[2] - m_min[2];
		FT length = (std::min(std::min(w, h), d)) / 5.0f;
		if (::glIsList(m_gl_bb_disp_list1)) {
			::glDeleteLists(m_gl_bb_disp_list1, 1);
		}
		m_gl_bb_disp_list1 = ::glGenLists(1);
		::glNewList(m_gl_bb_disp_list1, GL_COMPILE);
		::glPushAttrib(GL_ENABLE_BIT|GL_CURRENT_BIT);
		::glDisable(GL_LIGHTING);
		::glDisable(GL_TEXTURE_2D);
		::glColor3f(1.0f, 1.0f, 1.0f);
		::glBegin(GL_LINES);
		// Far Left Buttom
		::glVertex3f(m_min[0], m_min[1], m_min[2]);
		::glVertex3f(m_min[0]+length, m_min[1], m_min[2]);
		::glVertex3f(m_min[0], m_min[1], m_min[2]);
		::glVertex3f(m_min[0], m_min[1]+length, m_min[2]);
		::glVertex3f(m_min[0], m_min[1], m_min[2]);
		::glVertex3f(m_min[0], m_min[1], m_min[2]+length);
		// Far Right Buttom
		::glVertex3f(m_max[0], m_min[1], m_min[2]);
		::glVertex3f(m_max[0]-length, m_min[1], m_min[2]);
		::glVertex3f(m_max[0], m_min[1], m_min[2]);
		::glVertex3f(m_max[0], m_min[1]+length, m_min[2]);
		::glVertex3f(m_max[0], m_min[1], m_min[2]);
		::glVertex3f(m_max[0], m_min[1], m_min[2]+length);
		// Far Right Top
		::glVertex3f(m_max[0], m_max[1], m_min[2]);
		::glVertex3f(m_max[0]-length, m_max[1], m_min[2]);
		::glVertex3f(m_max[0], m_max[1], m_min[2]);
		::glVertex3f(m_max[0], m_max[1]-length, m_min[2]);
		::glVertex3f(m_max[0], m_max[1], m_min[2]);
		::glVertex3f(m_max[0], m_max[1], m_min[2]+length);
		// Far Left Top
		::glVertex3f(m_min[0], m_max[1], m_min[2]);
		::glVertex3f(m_min[0]+length, m_max[1], m_min[2]);
		::glVertex3f(m_min[0], m_max[1], m_min[2]);
		::glVertex3f(m_min[0], m_max[1]-length, m_min[2]);
		::glVertex3f(m_min[0], m_max[1], m_min[2]);
		::glVertex3f(m_min[0], m_max[1], m_min[2]+length);
		// Near Left Buttom
		::glVertex3f(m_min[0], m_min[1], m_max[2]);
		::glVertex3f(m_min[0]+length, m_min[1], m_max[2]);
		::glVertex3f(m_min[0], m_min[1], m_max[2]);
		::glVertex3f(m_min[0], m_min[1]+length, m_max[2]);
		::glVertex3f(m_min[0], m_min[1], m_max[2]);
		::glVertex3f(m_min[0], m_min[1], m_max[2]-length);
		// Near Right Buttom
		::glVertex3f(m_max[0], m_min[1], m_max[2]);
		::glVertex3f(m_max[0]-length, m_min[1], m_max[2]);
		::glVertex3f(m_max[0], m_min[1], m_max[2]);
		::glVertex3f(m_max[0], m_min[1]+length, m_max[2]);
		::glVertex3f(m_max[0], m_min[1], m_max[2]);
		::glVertex3f(m_max[0], m_min[1], m_max[2]-length);
		// Near Right Top
		::glVertex3f(m_max[0], m_max[1], m_max[2]);
		::glVertex3f(m_max[0]-length, m_max[1], m_max[2]);
		::glVertex3f(m_max[0], m_max[1], m_max[2]);
		::glVertex3f(m_max[0], m_max[1]-length, m_max[2]);
		::glVertex3f(m_max[0], m_max[1], m_max[2]);
		::glVertex3f(m_max[0], m_max[1], m_max[2]-length);
		// Near Left Top
		::glVertex3f(m_min[0], m_max[1], m_max[2]);
		::glVertex3f(m_min[0]+length, m_max[1], m_max[2]);
		::glVertex3f(m_min[0], m_max[1], m_max[2]);
		::glVertex3f(m_min[0], m_max[1]-length, m_max[2]);
		::glVertex3f(m_min[0], m_max[1], m_max[2]);
		::glVertex3f(m_min[0], m_max[1], m_max[2]-length);
		::glEnd();
		::glPopAttrib();
		::glEndList();
	}
	void gl_draw_bounding_box1()
	{
		glCallList(m_gl_bb_disp_list1);
	}

	void gl_build_bounding_box_disp_list2(void)  //Draw BBOX only at 8 corners
	{
		if (size_of_vertices() == 0) {
			return;
		}
		FT w, h, d;
		// calculate model width, height, and depth
		w = m_max[0] - m_min[0];
		h = m_max[1] - m_min[1];
		d = m_max[2] - m_min[2];
		FT length = (std::min(std::min(w, h), d)) / 7.0f;
		if (::glIsList(m_gl_bb_disp_list2)) {
			::glDeleteLists(m_gl_bb_disp_list2, 1);
		}
		m_gl_bb_disp_list2 = ::glGenLists(1);
		::glNewList(m_gl_bb_disp_list2, GL_COMPILE);
		::glPushAttrib(GL_ENABLE_BIT|GL_CURRENT_BIT);
		::glDisable(GL_LIGHTING);
		::glDisable(GL_TEXTURE_2D);
		::glColor3f(0.75f, 0.75f, 0.75f);
		::glDisable(GL_LINE_STIPPLE);
		::glLineWidth(4);
		::glBegin(GL_LINES);
		{
			// Far Left Buttom
			::glVertex3f(m_min[0], m_min[1], m_min[2]);
			::glVertex3f(m_min[0]+length, m_min[1], m_min[2]);
			::glVertex3f(m_min[0], m_min[1], m_min[2]);
			::glVertex3f(m_min[0], m_min[1]+length, m_min[2]);
			::glVertex3f(m_min[0], m_min[1], m_min[2]);
			::glVertex3f(m_min[0], m_min[1], m_min[2]+length);
			// Far Right Buttom
			::glVertex3f(m_max[0], m_min[1], m_min[2]);
			::glVertex3f(m_max[0]-length, m_min[1], m_min[2]);
			::glVertex3f(m_max[0], m_min[1], m_min[2]);
			::glVertex3f(m_max[0], m_min[1]+length, m_min[2]);
			::glVertex3f(m_max[0], m_min[1], m_min[2]);
			::glVertex3f(m_max[0], m_min[1], m_min[2]+length);
			// Far Right Top
			::glVertex3f(m_max[0], m_max[1], m_min[2]);
			::glVertex3f(m_max[0]-length, m_max[1], m_min[2]);
			::glVertex3f(m_max[0], m_max[1], m_min[2]);
			::glVertex3f(m_max[0], m_max[1]-length, m_min[2]);
			::glVertex3f(m_max[0], m_max[1], m_min[2]);
			::glVertex3f(m_max[0], m_max[1], m_min[2]+length);
			// Far Left Top
			::glVertex3f(m_min[0], m_max[1], m_min[2]);
			::glVertex3f(m_min[0]+length, m_max[1], m_min[2]);
			::glVertex3f(m_min[0], m_max[1], m_min[2]);
			::glVertex3f(m_min[0], m_max[1]-length, m_min[2]);
			::glVertex3f(m_min[0], m_max[1], m_min[2]);
			::glVertex3f(m_min[0], m_max[1], m_min[2]+length);
			// Near Left Buttom
			::glVertex3f(m_min[0], m_min[1], m_max[2]);
			::glVertex3f(m_min[0]+length, m_min[1], m_max[2]);
			::glVertex3f(m_min[0], m_min[1], m_max[2]);
			::glVertex3f(m_min[0], m_min[1]+length, m_max[2]);
			::glVertex3f(m_min[0], m_min[1], m_max[2]);
			::glVertex3f(m_min[0], m_min[1], m_max[2]-length);
			// Near Right Buttom
			::glVertex3f(m_max[0], m_min[1], m_max[2]);
			::glVertex3f(m_max[0]-length, m_min[1], m_max[2]);
			::glVertex3f(m_max[0], m_min[1], m_max[2]);
			::glVertex3f(m_max[0], m_min[1]+length, m_max[2]);
			::glVertex3f(m_max[0], m_min[1], m_max[2]);
			::glVertex3f(m_max[0], m_min[1], m_max[2]-length);
			// Near Right Top
			::glVertex3f(m_max[0], m_max[1], m_max[2]);
			::glVertex3f(m_max[0]-length, m_max[1], m_max[2]);
			::glVertex3f(m_max[0], m_max[1], m_max[2]);
			::glVertex3f(m_max[0], m_max[1]-length, m_max[2]);
			::glVertex3f(m_max[0], m_max[1], m_max[2]);
			::glVertex3f(m_max[0], m_max[1], m_max[2]-length);
			// Near Left Top
			::glVertex3f(m_min[0], m_max[1], m_max[2]);
			::glVertex3f(m_min[0]+length, m_max[1], m_max[2]);
			::glVertex3f(m_min[0], m_max[1], m_max[2]);
			::glVertex3f(m_min[0], m_max[1]-length, m_max[2]);
			::glVertex3f(m_min[0], m_max[1], m_max[2]);
			::glVertex3f(m_min[0], m_max[1], m_max[2]-length);
		}
		::glEnd();
		//middle segment using stipple style
		::glLineWidth(1);
		::glEnable(GL_LINE_STIPPLE);
		::glLineStipple(2,0x1C47);
		::glBegin(GL_LINES);
		{
			//4 segs parallel to X-axis
			::glVertex3f(m_min[0]+length, m_min[1], m_min[2]);
			::glVertex3f(m_max[0]-length, m_min[1], m_min[2]);
			::glVertex3f(m_min[0]+length, m_max[1], m_min[2]);
			::glVertex3f(m_max[0]-length, m_max[1], m_min[2]);
			::glVertex3f(m_min[0]+length, m_min[1], m_max[2]);
			::glVertex3f(m_max[0]-length, m_min[1], m_max[2]);
			::glVertex3f(m_min[0]+length, m_max[1], m_max[2]);
			::glVertex3f(m_max[0]-length, m_max[1], m_max[2]);
			//4 segs parallel to Y-axis
			::glVertex3f(m_min[0], m_min[1]+length, m_min[2]);
			::glVertex3f(m_min[0], m_max[1]-length, m_min[2]);
			::glVertex3f(m_min[0], m_min[1]+length, m_max[2]);
			::glVertex3f(m_min[0], m_max[1]-length, m_max[2]);
			::glVertex3f(m_max[0], m_min[1]+length, m_min[2]);
			::glVertex3f(m_max[0], m_max[1]-length, m_min[2]);
			::glVertex3f(m_max[0], m_min[1]+length, m_max[2]);
			::glVertex3f(m_max[0], m_max[1]-length, m_max[2]);
			//4 segs parallel to Z-axis
			::glVertex3f(m_min[0], m_min[1], m_min[2]+length);
			::glVertex3f(m_min[0], m_min[1], m_max[2]-length);
			::glVertex3f(m_min[0], m_max[1], m_min[2]+length);
			::glVertex3f(m_min[0], m_max[1], m_max[2]-length);
			::glVertex3f(m_max[0], m_min[1], m_min[2]+length);
			::glVertex3f(m_max[0], m_min[1], m_max[2]-length);
			::glVertex3f(m_max[0], m_max[1], m_min[2]+length);
			::glVertex3f(m_max[0], m_max[1], m_max[2]-length);
		}
		::glEnd();
		::glDisable(GL_LINE_STIPPLE);
		::glPopAttrib();
		::glEndList();
	}

	void gl_draw_bounding_box2(void)
	{
		::glCallList(m_gl_bb_disp_list2);
	}

	// copy bounding box
	void copy_bounding_box(Enriched_polyhedron<kernel,items> *pMesh)
	{
		m_min[0] = pMesh->xmin(); m_max[0] = pMesh->xmax();
		m_min[1] = pMesh->ymin(); m_max[1] = pMesh->ymax();
		m_min[2] = pMesh->zmin(); m_max[2] = pMesh->zmax();
	}

	// degree of a face: #edges of this face, 3 for triangular mesh
	static unsigned int degree(Facet_handle pFace)
	{
		return CGAL::circulator_size(pFace->facet_begin());    
	}

	// valence of a vertex: # edges share this vertex
	static unsigned int valence(Vertex_handle pVertex)
	{
		return CGAL::circulator_size(pVertex->vertex_begin());
	}

	unsigned int fvid(Facet_handle pFace, int v)
	{
		Halfedge_around_facet_circulator pHE;
		pHE = pFace->facet_begin();
		if (v==0) {
			return pHE->vertex()->id();
		}
		pHE++;
		if (v==1) {
			return pHE->vertex()->id();
		}
		pHE++;
		if (v==2) {
			return pHE->vertex()->id();
		}
		return -1;
	}

	unsigned int fvid(int f, int v)
	{
		Halfedge_around_facet_circulator pHE;
		pHE = FH(f)->facet_begin();
		if (v==0) {
			return pHE->vertex()->id();
		}
		pHE++;
		if (v==1) {
			return pHE->vertex()->id();
		}
		pHE++;
		if (v==2) {
			return pHE->vertex()->id();
		}
		return -1;
	}

	const Point_3& fvp(Facet_handle pFace, int v)
	{
		Halfedge_around_facet_circulator pHE;
		pHE = pFace->facet_begin();
		if (v==0) {
			return pHE->vertex()->point();
		}
		pHE++;
		if (v==1) {
			return pHE->vertex()->point();
		}
		pHE++;
		return pHE->vertex()->point();
	}

	const Point_3& fvp(int f, int v)
	{
		Halfedge_around_facet_circulator pHE;
		pHE = FH(f)->facet_begin();
		if (v==0) {
			return pHE->vertex()->point();
		}
		pHE++;
		if (v==1) {
			return pHE->vertex()->point();
		}
		pHE++;
		return pHE->vertex()->point();
	}

	// get the id of the opposite facet of v (0, 1 or 2) in facet f
	int fvofid(int f, int v)
	{
		Halfedge_around_facet_circulator pHE;
		pHE = FH(f)->facet_begin();
		if (v==1) {
			if (pHE->opposite()->is_border()) {
				return -1;
			} else {
				return pHE->opposite()->facet()->id();
			}
		} else if (v==2) {
			if (pHE->next()->opposite()->is_border()) {
				return -1;
			} else {
				return pHE->next()->opposite()->facet()->id();
			}
		} else if (v==0) {
			if (pHE->next()->next()->opposite()->is_border()) {
				return -1;
			} else {
				return pHE->next()->next()->opposite()->facet()->id();
			}
		} else {
			return -1;
		}
	}

	// check whether	a vertex is on a boundary or not
	static bool	is_border(Vertex_handle	pVertex)
	{
		Halfedge_around_vertex_circulator	pHalfEdge	=	pVertex->vertex_begin();
		if(pHalfEdge ==	NULL)	// isolated	vertex
			return true;
		Halfedge_around_vertex_circulator	d	=	pHalfEdge;
		CGAL_For_all(pHalfEdge,d)
			if(pHalfEdge->is_border())
				return true;
		return false;
	}

	// get any border halfedge attached	to a vertex
	Halfedge_handle	get_border_halfedge(Vertex_handle pVertex)
	{
		Halfedge_around_vertex_circulator pHalfEdge = pVertex->vertex_begin();
		Halfedge_around_vertex_circulator d = pHalfEdge;
		CGAL_For_all(pHalfEdge,d)
			if(pHalfEdge->is_border())
				return pHalfEdge;
		return NULL;
	}

	// tag all halfedges
	void tag_halfedges(const int tag)
	{
		for(Halfedge_iterator pHalfedge = halfedges_begin();
			pHalfedge != halfedges_end();
			pHalfedge++)
			pHalfedge->tag(tag);
	}

	// tag all facets
	void tag_facets(const int tag)
	{
		for(Facet_iterator pFace=facets_begin(); pFace!=facets_end(); pFace++) {
			pFace->tag(tag);
		}
	}

	// tag all vertices
	void tag_vertices(const int tag)
	{
		for (Vertex_iterator pVertex=vertices_begin(); pVertex!=vertices_end(); pVertex++) {
			pVertex->tag(tag);
		}
	}

	// set index for all vertices: vertex->index / index->vertex
	void set_id_to_items_map()
	{
		index_to_vertex_map = Random_access_vertex_index(vertices_begin(), vertices_end());
		index_to_facet_map = Random_access_facet_index(facets_begin(), facets_end());
		index_to_halfedge_map = Random_access_halfedge_index(halfedges_begin(), halfedges_end());
		m_is_set_id_to_items_map = true;
	}

	void set_hds_items_id()
	{
		int	index = 0;
		for (Vertex_iterator pVertex=vertices_begin(); pVertex!=vertices_end(); pVertex++) {
			pVertex->id() = index++;
		}
		index = 0;
		for (Halfedge_iterator pHalfedge=halfedges_begin(); pHalfedge!=halfedges_end(); pHalfedge++) {
			pHalfedge->id() = index++;
		}
		index = 0;
		for (Face_iterator pFacet=facets_begin(); pFacet!=facets_end(); pFacet++) {
			pFacet->id() = index++;
		}
		m_is_set_items_id = true;
	}

	// is pure degree ?
	bool is_pure_degree(unsigned int d)
	{
		for(Facet_iterator pFace	=	facets_begin();
			pFace	!= facets_end();
			pFace++)
			if(degree(pFace) != d)
				return false;
		return true;
	}

	// compute type
	void compute_type()
	{
		m_pure_quad = is_pure_degree(4);
		m_pure_triangle = is_pure_degree(3);
	}

	// compute facet center
	void compute_facet_center(Facet_handle pFace, Point& center)
	{
		Halfedge_around_facet_circulator pHalfEdge = pFace->facet_begin();
		Halfedge_around_facet_circulator end = pHalfEdge;
		Vector vec(0.0,0.0,0.0);
		int	degree = 0;
		CGAL_For_all(pHalfEdge,end)
		{
			vec	= vec + (pHalfEdge->vertex()->point()-CGAL::ORIGIN);
			degree++;
		}
		center = CGAL::ORIGIN + (vec/(kernel::FT)degree);
	}

	// compute average edge length around a vertex
	FT average_edge_length_around(Vertex_handle pVertex)
	{
		FT sum = 0.0;
		Halfedge_around_vertex_circulator pHalfEdge = pVertex->vertex_begin();
		Halfedge_around_vertex_circulator end = pHalfEdge;
		Vector vec(0.0,0.0,0.0);
		int	degree = 0;
		CGAL_For_all(pHalfEdge,end)
		{
			Vector vec = pHalfEdge->vertex()->point()-
				pHalfEdge->opposite()->vertex()->point();
			sum += std::sqrt(vec*vec);
			degree++;
		}
		return sum / (FT) degree;
	}

	// compute minimum edge length around a vertex
	FT min_edge_length_around(Vertex_handle pVertex)
	{
		FT min_edge_length = 1e38;
		Halfedge_around_vertex_circulator pHalfEdge = pVertex->vertex_begin();
		Halfedge_around_vertex_circulator end = pHalfEdge;
		Vector vec(0.0,0.0,0.0);
		int	degree = 0;
		CGAL_For_all(pHalfEdge,end)
		{
			Vector vec = pHalfEdge->vertex()->point()-
				pHalfEdge->opposite()->vertex()->point();
			FT len = std::sqrt(vec*vec);
			if(len < min_edge_length)
				min_edge_length = len;
		}
		return min_edge_length;
	}

	bool pick_vertex_by_idx(unsigned vid)
	{
		if(vid<0 || vid> size_of_vertices()-1)
			return false;
		Vertex_iterator pVertex = VH(vid);
		if(!pVertex->have_tag(PICKED)){
			picked_vertices.push_back(vid);
			pVertex->add_tag(PICKED);
			return true;
		}
		else	return false;		
	}

	bool pick_facet_by_idx(unsigned fid)
	{
		if(fid<0 || fid > size_of_facets()-1)
			return false;
		Facet_iterator pFacet = FH(fid);
		if(!pFacet->have_tag(PICKED)){
			picked_facets.push_back(fid);
			pFacet->add_tag(PICKED);
			return true;
		}
		else	return false;		
	}

	bool pick_edge_by_idx(unsigned eid)
	{
		if(eid<0 || eid > size_of_halfedges()-1)
			return false;
		Halfedge_iterator pHE = EH(eid);
		if(!pHE->have_tag(PICKED)){
			picked_edges.push_back(eid);
			pHE->add_tag(PICKED);
			return true;
		}
		else	return false;		
	}

	int pick_vertex_by_ray(Vector3 &sp, Vector3 &d, float epsilon,bool bMuiltipick=false,bool bUnpick = false)
	{
		int pi = -1;
		float min_t = FLT_MAX;
		float min_dist = epsilon, dist;
		float t = 0.0f;
		Vector3 vert, vec;

		d.normalize();
		Vertex_iterator pVertex;
		for(pVertex=vertices_begin(); pVertex!=vertices_end(); pVertex++)
		{	
			vert.set(pVertex->point()[0], pVertex->point()[1], pVertex->point()[2]);
			t = d.dot(vert - sp);
			if (t <= 0 || t >= min_t) {
				continue;
			}
			vec = vert - (sp + d * t);
			dist = vec.sqauremagnitude();
			if (dist < min_dist) {
				pi = pVertex->id();	// vertex index
				min_dist = dist;
				min_t = t;
			}
		}
		if (pi != -1) {	// picked
			pVertex = index_to_vertex_map[pi];
			if (!bUnpick) {
				if(!pVertex->have_tag(PICKED)){
					picked_vertices.push_back(pi);
					pVertex->add_tag(PICKED);	
				}				
			} else {
				if(pVertex->have_tag(PICKED)){
					picked_vertices.remove(pi);
					pVertex->remove_tag(PICKED);
				}
			}
			return pi;
		}
		return -1;
	}

	int pick_edge_by_ray(Vector3 &sp, Vector3 &d, float epsilon,bool bMuiltipick=false,bool bUnpick = false)
	{
		int Eidx = -1;  
		float min_t = FLT_MAX;
		float min_dist = epsilon, dist;
		float t = 0.0f;
		Vector3 vec;

		d.normalize();
		Vector3 E_mid_point; //Vector E_mid_point_p;
		Halfedge_iterator pE = halfedges_begin();
		CGAL_For_all(pE, halfedges_end())
		{
			const Point& p1 = pE->vertex()->point();
			const Point& p2 = pE->opposite()->vertex()->point();
			E_mid_point.set((p1[0] + p2[0])*0.5,(p1[1] + p2[1])*0.5,(p1[2] + p2[2])*0.5);
			t = d.dot(E_mid_point - sp);
			if (t <= 0 || t >= min_t) 	
				continue;
			vec = E_mid_point - (sp + d * t);
			dist = vec.sqauremagnitude();
			if (dist < min_dist) {
				Eidx = pE->id();	// edge index
				min_dist = dist;
				min_t = t;
			}
		}
		if (Eidx != -1) {	// picked
			pE = index_to_halfedge_map[Eidx];
			if (!bUnpick) {
				if(!pE->have_tag(PICKED)){
					picked_edges.push_back(Eidx);
					pE->add_tag(PICKED);
				}
			} else {
				picked_edges.remove(Eidx);
				pE->remove_tag(PICKED);
			}
			return Eidx;
		}
		return -1;
	}

	void clear_picked_edges(void)
	{
		Halfedge_iterator eh;
		for (std::list<int>::iterator it =picked_edges.begin(); it!=picked_edges.end(); ++it){
			eh = index_to_halfedge_map[(*it)];
			eh->remove_tag(PICKED);
		}
		picked_edges.clear();
	}

	int nb_picked_edges(void)
	{
		return picked_edges.size();
	}

	int pick_vertex_by_ray_normal(Vector3 &sp, Vector3 &d, float epsilon)
	{
		int pi = -1;
		float min_t = FLT_MAX;
		float min_dist = epsilon, dist;
		float t = 0.0f;
		Vector3 vert;
		d.normalize();
		for(Vertex_iterator pVertex=vertices_begin(); pVertex!=vertices_end(); pVertex++)
		{
			vert.set(pVertex->point()[0], pVertex->point()[1], pVertex->point()[2]);
			t = d.dot(vert - sp);
			if (t <= 0 || t >= min_t) {
				continue;
			}
			vert = vert - (sp + d * t);
			dist = vert.sqauremagnitude();
			if (dist < min_dist) {
				vert.set(pVertex->normal()[0], pVertex->normal()[1], pVertex->normal()[2]);
				if (vert.dot(d) < 0.2) {
					pi = pVertex->id();	// vertex index
					min_dist = dist;
					min_t = t;
				}
			}
		}
		return pi;
	}

	void pick_mh_hdl_by_ray_n(const Vector3 &sp, const Vector3 &d, int &ph, int &pv, float epsilon)
	{
		float min_t=FLT_MAX;
		float min_dist=epsilon, dist;
		float t=0.0f;
		Vector3 vert;
		Vector3 dir(d);
		dir.normalize();
		std::set<int>::iterator v;
		ph = -1;
		pv = -1;
		for (int i=0; i<m_hdl_v.size(); i++)
		{
			for(v=m_hdl_v[i].begin(); v!=m_hdl_v[i].end(); v++)
			{
				vert.set(VP(*v)[0], VP(*v)[1], VP(*v)[2]);
				t = dir.dot(vert - sp);
				if (t <= 0 || t >= min_t) {
					continue;
				}
				vert = vert - (sp + dir * t);
				dist = vert.sqauremagnitude();
				if (dist < min_dist) {
					vert.set(VN(*v)[0], VN(*v)[1], VN(*v)[2]);
					if (vert.dot(dir) < 0) {
						ph = i;
						pv = *v;	// vertex index
						min_dist = dist;
						min_t = t;
					}
				}
			}
		}
	}

	int pick_vertices_by_column(Vector3 &sp, Vector3 &dx, Vector3 &dy, Vector3 &dz, bool bFrontOnly=false, bool bSilhouetteOnly=false, FT viewOscLow=-0.0, FT viewOscHigh=1000.0, FT viewFrontLow=0.0)
	{
		float t1 = 0.0f, t2 = 0.0f, t3 = 0.0f;

		float width = dx.magnitude();
		float height = dy.magnitude();
		dx.normalize();
		dy.normalize();
		dz.normalize();

		Vector3 vec, vert, norml;
		int last_found_Vidx = -1;
		for(Vertex_iterator pVertex=vertices_begin(); pVertex!=vertices_end(); pVertex++)
		{
			vert.set(pVertex->point()[0], pVertex->point()[1], pVertex->point()[2]);
			vec = vert - sp;
			if (dx.dot(vec)<=0 || dy.dot(vec)<=0 || dz.dot(vec)<=0) {
				continue;
			}
			if (width < fabs(dx.dot(vec))) {
				continue;
			}
			if (height < fabs(dy.dot(vec))) {
				continue;
			}
			if(bFrontOnly){
				norml.set(pVertex->normal()[0], pVertex->normal()[1], pVertex->normal()[2]);
				if (norml.dot(vec)>viewFrontLow)	continue;
			}
			else if(bSilhouetteOnly){
				norml.set(pVertex->normal()[0], pVertex->normal()[1], pVertex->normal()[2]);
				if (norml.dot(vec)<=viewOscLow || norml.dot(vec)>=viewOscHigh)	continue;
			}
			last_found_Vidx = pVertex->id();
			if(pVertex->have_tag(PICKED)){ //already picked
				continue;
			}
			pVertex->add_tag(PICKED);
			picked_vertices.push_back(pVertex->id());
		}
		return last_found_Vidx;
	}

	void unpick_vertices_by_column(Vector3 &sp, Vector3 &dx, Vector3 &dy, Vector3 &dz, bool bFrontOnly=false, bool bSilhouetteOnly=false, FT viewOscLow=-0.0, FT viewOscHigh=1000.0, FT viewFrontLow=0.0)
	{
		float t1 = 0.0f, t2 = 0.0f, t3 = 0.0f;
		float width = dx.magnitude();
		float height = dy.magnitude();
		dx.normalize();
		dy.normalize();
		dz.normalize();
		Vector3 vec, vert, norml;
		std::vector<std::list<int>::iterator> its;
		Vertex_iterator pV;
		for (std::list<int>::iterator it=picked_vertices.begin(); it!=picked_vertices.end(); it++)
		{
			pV = index_to_vertex_map[(*it)];
			vert.set(pV->point().x(), pV->point().y(), pV->point().z());
			vec = vert - sp;
			if (dx.dot(vec)<=0 || dy.dot(vec)<=0 || dz.dot(vec)<=0) {
				continue;
			}
			if (width < fabs(dx.dot(vec))) {
				continue;
			}
			if (height < fabs(dy.dot(vec))) {
				continue;
			}
			if(bFrontOnly){
				norml.set(pV->normal()[0], pV->normal()[1], pV->normal()[2]);
				if (norml.dot(vec)>viewFrontLow)	continue;
			}
			else if(bSilhouetteOnly){
				norml.set(pV->normal()[0], pV->normal()[1], pV->normal()[2]);
				if (norml.dot(vec)<=viewOscLow || norml.dot(vec)>=viewOscHigh)	continue;
			}
			its.push_back(it);
			pV->remove_tag(PICKED);
		}
		for (unsigned int i=0; i<its.size(); i++)
		{
			picked_vertices.erase(its[i]);
		}
	}

	void clear_picked_vertices(void)
	{
		Vertex_iterator vh;
		for (std::list<int>::iterator it =picked_vertices.begin(); it!=picked_vertices.end(); ++it){
			vh = index_to_vertex_map[(*it)];
			vh->remove_tag(PICKED);
		}
		picked_vertices.clear(); 
	}

	int nb_picked_vertices(void)
	{
		return picked_vertices.size();
	}
	void write_picked_vertices(std::ofstream &stream)
	{
		stream << picked_vertices.size() << std::endl;
		std::set<int>::iterator pv;
		for (pv=picked_vertices.begin(); pv!=picked_vertices.end(); pv++)
		{
			stream << *pv << std::endl;
		}
	}

	void write_picked_facets(std::ofstream &stream)
	{
		stream << picked_facets.size() << std::endl;
		std::list<int>::iterator pf;
		for (pf=picked_facets.begin(); pf!=picked_facets.end(); pf++)
		{
			stream << *pf << std::endl;
		}
	}

	void read_picked_vertices(std::ifstream &stream)
	{
		int nV, v;
		stream >> nV;
		int nMV = size_of_vertices();
		if (nMV <= nV)
		{
			return;
		}
		picked_vertices.clear();
		for (int i=0; i<nV; i++) {
			stream >> v;
			if (v >= nMV)
			{
				v = nMV-1;
			}
			picked_vertices.push_back(v);
		}
	}

	void read_picked_facets(std::ifstream &stream)
	{
		int nF, f;
		stream >> nF;
		int nMF = size_of_facets();
		if (nMF <= nF)
		{
			return;
		}
		picked_facets.clear();
		for (int i=0; i<nF; i++) {
			stream >> f;
			if (f >= nMF)
			{
				f = nMF-1;
			}
			picked_facets.push_back(f);
		}
	}
	void Write_Texcoords(std::ofstream &stream)
	{
		int	nV;
		stream << "Texcoords" << std::endl; 
		nV= size_of_vertices();		stream << nV << std::endl; 
		if(nV!=0){
			Vertex_const_iterator it = vertices_begin();
			Vertex_const_iterator it_end = vertices_end();
			CGAL_For_all(it , it_end){
				stream<<it->u()<<","<<it->v()<<std::endl;
			}			
		}
	}
	void drag_picked_vertices(CPoint move, Matrix3 rot_matx,GLdouble *modelview, GLdouble *projection, GLint *viewport)
	{
		std::list<int>::const_iterator it =picked_vertices.begin();
		Vector3 move3d = rot_matx*Vector3(move.x,-move.y,0.0);
		Vector mov(move3d.x*0.003,move3d.y*0.003,move3d.z*0.003);
		for (; it!=picked_vertices.end(); ++it){
			VP(*it) = VP(*it) + mov; 	
		}
	}
	bool Load_Texcoords(std::ifstream &stream)
	{
		int l(0);
		int	nV;
		std::string buf;
		char s;
		stream >> buf; l++;
		if(buf!="Texcoords"){
			MSG_BOX_ERROR("load texcoords : error type head at %d line.", l);		return false;
		}
		stream>>nV; l++;
		if(nV != size_of_vertices() ){
			MSG_BOX_ERROR("load texcoords : vertices' number %d not equal to size of vertices at %d line.", nV);		
			return false;
		}
		FT u,v;
		Vertex_iterator it = vertices_begin();	
		Vertex_iterator it_end = vertices_end();
		CGAL_For_all(it , it_end){
			stream>>u>>s>>v;
			it->texcoord(u,v);
		}			
		return true;		
	}

	void Copy_Points()
	{
		if(!m_bPointsCopied){
			m_points_copy.clear();	m_points_copy.resize(size_of_vertices());
			for (unsigned i=0; i<size_of_vertices(); ++i)
				m_points_copy[i].set(VP(i).x(),VP(i).y(),VP(i).z());
			m_bPointsCopied =true;
		}
	}
	void CopyBack_Points()
	{
		if(m_bPointsCopied){
			for (unsigned i=0; i<size_of_vertices(); ++i)
				VP(i)  = Point(m_points_copy[i].x,m_points_copy[i].y,m_points_copy[i].z);
		}
	}

	void Rotate_model(Quaternion<> &rot)
	{	
		if(m_bPointsCopied){
			Vector3 new_p;
			for (unsigned i=0; i<size_of_vertices(); ++i){
				Point_3 &p=VP(i);
				new_p = rot*m_points_copy[i];
				p = Point_3(new_p.x,new_p.y,new_p.z);
			}
		}
	}
	void vertices_2d_projection(std::vector<Vector2> &v2d,GLdouble *modelview, GLdouble *projection, GLint *viewport)
	{
		int	nV;	nV= size_of_vertices(); 
		if(v2d.size()!=nV){
			v2d.clear();v2d.resize(nV);
		}		
		/*CAUTION:  if any doubt about how Opengl manipulate the modelview matrix,.. 		
		std::cout<<"\nroation:  \n"; 
		for (int i=0; i<3; ++i)		for(int j=0; j<3; ++j)
		std::cout<<modelview[i*4+j]/rotation.M[i*3+j]<<" ";
		std::cout<<"\ntranslation:  \n"; 
		for (int i=12; i<16; ++i)	std::cout<<modelview[i]<<" ";
		*/			
		double winz;
		for (int vi=0; vi<nV; ++vi){  
			::gluProject(VP(vi)[0], VP(vi)[1], VP(vi)[2],
				modelview, projection, viewport, &v2d[vi].x, &v2d[vi].y, &winz);
			/*if you want the coord system origin point at the middle....*/
			//p_2d.x -= (window_size.cx*0.5); p_2d.y -= (window_size.cy*0.5);				
		}
	}
	void vnormals_2d_projection(std::vector<Vector2> &vn2d, Matrix3 &view_rot)
	{
		int	nV;	nV= size_of_vertices(); 
		if(vn2d.size()!=nV){
			vn2d.clear();vn2d.resize(nV);
		}
		//Matrix3	inv_rotation = inv_rotation.invert(rotation);
		view_rot.transpose(); //it's the same as .invert for the orthogonal matrix.
		unsigned i=0;
		Vector3 VN3d;
		for (int vi=0; vi<nV; ++vi){
			VN3d = view_rot * Vector3(VN(vi)[0], VN(vi)[1], VN(vi)[2]);		
			vn2d[i] = Vector2(VN3d.x, VN3d.y);
			i++;
		}
	}
	void picked_vertices_2d_projection(std::vector<Vector2> &v2d, GLdouble *modelview, GLdouble *projection, GLint *viewport)
	{
		int	nPickedV;	nPickedV= nb_picked_vertices(); 
		if(v2d.size() != nPickedV){
			v2d.clear();v2d.resize(nPickedV);
		}		
		double winz;
		std::list<int>::const_iterator it = picked_vertices.begin();
		int i=0;
		for (; it!=picked_vertices.end(); ++it)
		{
			::gluProject(VP(*it)[0], VP(*it)[1], VP(*it)[2],
				modelview, projection, viewport, &v2d[i].x, &v2d[i].y, &winz);
			++i;
		}
	}
	void picked_vnormals_2d_projection(std::vector<Vector2> &vn2d, Matrix3 &view_rot)
	{
		int	nV;	nV= nb_picked_vertices(); 
		if(vn2d.size()!=nV){
			vn2d.clear();vn2d.resize(nV);
		}
		std::list<int>::const_iterator it = picked_vertices.begin();
		view_rot.transpose();
		unsigned i=0;
		Vector3 VN3d;
		for (; it!=picked_vertices.end(); ++it)	{
			VN3d = view_rot * Vector3(VN(*it)[0], VN(*it)[1], VN(*it)[2]);		
			vn2d[i] = Vector2(VN3d.x, VN3d.y);
			i++;
		}
	}
	void Save_Vertices_and_Normal_2DProjection(std::ofstream &stream, Matrix3 &rotation, GLdouble *modelview, GLdouble *projection, GLint *viewport)
	{
		int	nV;
		stream << "Vertex Projection" << std::endl; 
		nV= size_of_vertices();		stream << nV << std::endl; 

		rotation.transpose();

		Vector3 p_2d(.0,.0,.0), VN_2d(.0,.0,.0);
		// all mesh V and VNormal's screen version.
		for (int vi=0; vi<nV; ++vi){  
			::gluProject(VP(vi)[0], VP(vi)[1], VP(vi)[2],
				modelview, projection, viewport,
				&p_2d.x, &p_2d.y, &p_2d.z);
			/*if you want the coord system origin point at the middle....*/
			//p_2d.x -= (window_size.cx*0.5); p_2d.y -= (window_size.cy*0.5);
			VN_2d = rotation * Vector3(VN(vi)[0], VN(vi)[1], VN(vi)[2]);	
			stream << p_2d << VN_2d << std::endl;			
		}
	}

	void Write_Pick_Objects(std::ofstream &stream)
	{
		int	nV, nF, nE;	
		std::list<int>::const_iterator it;

		stream << "Vertex" << std::endl; 
		nV= nb_picked_vertices();		stream << nV << std::endl; 
		if(nV!=0){
			it = picked_vertices.begin();
			for (;it!=picked_vertices.end(); ++it)
				stream<<(*it)<<" ";
			stream <<std::endl;
		}

		stream << "Face" << std::endl; 
		nF= nb_picked_facets();		stream << nF << std::endl; 
		if(nF!=0){
			it = picked_facets.begin();
			for (;it!=picked_facets.end(); ++it)
				stream<<(*it)<<" ";
			stream <<std::endl;
		}

		stream << "Edge" << std::endl; 
		nE= nb_picked_edges();		stream << nE << std::endl; 
		if(nE!=0){
			it = picked_edges.begin();
			for (;it!=picked_edges.end(); ++it)
				stream<<(*it)<<" ";
		}
	}

	void Load_Pick_Objects(std::ifstream &stream)
	{
		int l(0);
		int	nV, nF,nE;
		int v,f,e;
		std::string buf;

		stream >> buf; l++;
		if(buf!="Vertex"){
			MSG_BOX_ERROR("load picked objects : error vertex part head at %d line.", l);		return;
		}
		stream>>nV; l++;
		if(nV != 0 ){
			if( !ValidRange(nV,0,size_of_vertices()) ){
				MSG_BOX_ERROR("load picked objects : vertices' number not valid range at %d line.", l);		return;
			}
			picked_vertices.clear();
			tag_vertices(0);
			for (int i=0; i<nV; ++i){
				stream>>v; l++;
				picked_vertices.push_back(v);
			}
		}

		stream >> buf; l++;
		if(buf!="Face"){
			MSG_BOX_ERROR("load picked objects : error facet part head at %d line.", l);		return;
		}
		stream>>nF; l++;
		if(nF != 0 ){
			if( !ValidRange(nF,0,size_of_facets()) ){
				MSG_BOX_ERROR("load picked objects : facets' number not valid range at %d line.", l);		return;
			}
			picked_facets.clear();
			tag_facets(0);
			for (int i=0; i<nF; ++i){
				stream>>f; l++;			
				picked_facets.push_back(f);
			}
		}

		stream >> buf; l++;
		if(buf!="Edge"){
			MSG_BOX_ERROR("load picked objects : error edge part head at %d line.", l);		return;
		}
		stream>>nE; l++;
		if(nE!=0){
			if( !ValidRange(nE,0,size_of_halfedges()) ){
				MSG_BOX_ERROR("load picked objects : edges' number not valid range at %d line.", l);		return;
			}
			picked_edges.clear();
			tag_halfedges(0);
			for (int i=0; i<nE; ++i){
				stream>>e; 		
				picked_edges.push_back(e);
			}
		}
	}
	void suitable_surface_point_size(FT fGLWidth, CSize window_size, FT fCamZoom=1.0)
	{		
		FT dR = get_avg_edge_len();
		m_suitable_point_size = dR * 0.2 * window_size.cx/fGLWidth *fCamZoom;
	}

	void gl_draw_picked_vertices(FT fGLWidth, CSize window_size, FT fCamZoom=1.0)
	{
		if (picked_vertices.size() == 0) 	return;

		suitable_surface_point_size(fGLWidth,  window_size, fCamZoom);
		//::glPushMatrix();
		::glPushAttrib(GL_ALL_ATTRIB_BITS);
		::glDisable(GL_LIGHTING);
		::glDisable(GL_TEXTURE_2D);
		::glColor4f(1.0f, 0.0f, 0.0f,1.0);
		std::list<int>::iterator pv;
		Vertex_handle v;
		::glDisable(GL_BLEND);
		//::glBlendFunc(GL_ONE,GL_SRC_ALPHA_SATURATE);
		//glPointSize(m_suitable_point_size);
		//::glBegin(GL_POINTS);		
		for (pv=picked_vertices.begin(); pv!=picked_vertices.end(); pv++) {
			v = index_to_vertex_map[(*pv)];
			::glPushMatrix();
			::glTranslatef(v->halfedge()->vertex()->point().x(),
							v->halfedge()->vertex()->point().y(),
							v->halfedge()->vertex()->point().z());
			::glutSolidSphere(0.014f, 20, 20);
			::glPopMatrix();
// 			::glVertex3f(v->halfedge()->vertex()->point().x(), 
// 				v->halfedge()->vertex()->point().y(), 
// 				v->halfedge()->vertex()->point().z());
		}
		//::glEnd();
		//::glPopMatrix();
		::glPopAttrib();
	}

	void gl_draw_picked_edges(void)
	{
		if (picked_edges.size() == 0) {
			return;
		}
		std::list<int>::iterator it_pickededge;
		Halfedge_iterator pHE;
		::glPushAttrib(GL_ALL_ATTRIB_BITS/*GL_LIGHTING_BIT|GL_CURRENT_BIT*/);
		::glDisable(GL_LIGHTING);
		::glDisable(GL_TEXTURE_2D);
		GLfloat color[3] = {0.25f, 0.87f, 0.15f}; 
		glColor3fv(color);
		::glLineWidth(3.5);
		::glBegin(GL_LINES);	
		for (it_pickededge=picked_edges.begin(); it_pickededge!=picked_edges.end(); it_pickededge++) {
			pHE = index_to_halfedge_map[(*it_pickededge)];
			glVertex3f(pHE->vertex()->point().x(),pHE->vertex()->point().y(),pHE->vertex()->point().z());
			glVertex3f(pHE->opposite()->vertex()->point().x(),pHE->opposite()->vertex()->point().y(),pHE->opposite()->vertex()->point().z());
		}
		::glEnd();
		::glPopAttrib();
	}

	int intersect_with_ray(FT *startp, FT *dir)
	{
		Vector3 sp(startp), d(dir);
		Vector3 e1, e2, p, s, q;
		float   t(0.0f), u(0.0f), v(0.0f), w(0.0f), tmp(0.0f);
		Point_3 v1, v2, v3;
		int     times = 0;
		Halfedge_around_facet_circulator pHalfedge;

		for (Facet_iterator pFacet=facets_begin(); pFacet!=facets_end(); pFacet++)
		{
			pHalfedge = pFacet->facet_begin();
			v1 = pHalfedge->vertex()->point();
			v2 = pHalfedge->next()->vertex()->point();
			v3 = pHalfedge->next()->next()->vertex()->point();
			e1.set(v2.x()-v1.x(), v2.y()-v1.y(), v2.z()-v1.z());
			e2.set(v3.x()-v1.x(), v3.y()-v1.y(), v3.z()-v1.z());
			p = d.cross(e2);
			tmp = p.dot(e1);
			if (tmp > -FLT_EPSILON && tmp < FLT_EPSILON) {//0.000001f
				continue;
			}
			tmp = 1.0f / tmp;
			s.set(sp.x-v1.x(), sp.y-v1.y(), sp.z-v1.z());
			u = tmp * p.dot(s);
			if (u < 0.0f || u > 1.0f) {
				continue;
			}
			q = s.cross(e1);
			v = tmp * q.dot(d);
			if (v < 0.0f || v > 1.0f) {
				continue;
			}
			w = u + v;
			if (w > 1.0f) {
				continue;
			}
			t = tmp * q.dot(e2);
			if (t > 0) {
				times++;
			}
		}
		return times;
	}

	void check_picked(int f)
	{
		Facet_handle FH = index_to_facet_map[f];
		Halfedge_around_facet_circulator pHE = FH->facet_begin();
		int v1, v2, ov1, ov2;
		do {
			v1 = pHE->vertex()->id();
			v2 = pHE->next()->vertex()->id();
			ov1 = pHE->opposite()->next()->next()->vertex()->id();
			ov2 = pHE->opposite()->vertex()->id();
		} while(++pHE != FH->facet_begin());
	}
	int pick_facet_by_ray(Vector3 &sp, Vector3 &d,bool bMuiltipick=false, bool bUnpick = false)
	{
		Vector3 e1, e2, p, s, q;
		FT   t(0), u(0), v(0), w(0), tmp(0);
		FT   min_t = FLT_MAX;
		int     pi = -1;
		Point_3 v1, v2, v3;
		int		i = 0;
		Halfedge_around_facet_circulator pHalfedge;
		Facet_iterator pFacet;
		for(pFacet=facets_begin(); pFacet!=facets_end(); pFacet++, i++)
		{
			pHalfedge = pFacet->facet_begin();
			v1 = pHalfedge->vertex()->point();
			v2 = pHalfedge->next()->vertex()->point();
			v3 = pHalfedge->next()->next()->vertex()->point();
			e1.set(v2.x()-v1.x(), v2.y()-v1.y(), v2.z()-v1.z());
			e2.set(v3.x()-v1.x(), v3.y()-v1.y(), v3.z()-v1.z());
			p = d.cross(e2);
			tmp = p.dot(e1);

			if (tmp > -FLT_EPSILON && tmp < FLT_EPSILON) {
				continue;
			}
			tmp = 1.0f / tmp;
			s.set(sp.x-v1.x(), sp.y-v1.y(), sp.z-v1.z());
			u = tmp * p.dot(s);
			if (u < 0.0f || u > 1.0f) {
				continue;
			}
			q = s.cross(e1);
			v = tmp * q.dot(d);
			if (v < 0.0f || v > 1.0f) {
				continue;
			}
			w = u + v;
			if (w > 1.0f) {
				continue;
			}
			t = tmp * q.dot(e2);
			if (t > 0 && t < min_t) {
				pi = i;
				min_t = t;
			}
		}
		if (pi != -1) {
			pFacet = index_to_facet_map[pi];
			if (!bUnpick) {
				if(!pFacet->have_tag(PICKED)){
					picked_facets.push_back(pi);
					pFacet->add_tag(PICKED);	
				}				
			} else {
				if(pFacet->have_tag(PICKED)){
					picked_facets.remove(pi);
					pFacet->remove_tag(PICKED);
				}
			}
			return pi;
		}
		return -1;
	}

	int pick_facets_by_column(Vector3 &sp, Vector3 &dx, Vector3 &dy, Vector3 &dz, bool bFrontOnly=false, bool bSilhouetteOnly=false, FT viewOscLow=-0.0, FT viewOscHigh=1000.0, FT viewFrontLow=0.0)
	{
		FT t1(0), t2(0), t3(0);
		FT width = dx.magnitude();
		FT height = dy.magnitude();
		dx.normalize();
		dy.normalize();
		dz.normalize();
		Vector3 vec, norml;
		Halfedge_handle pHalfedge;
		Point_3 v1, v2, v3;
		int last_found_Fidx= 0;
		for(Facet_iterator pFacet=facets_begin(); pFacet!=facets_end(); pFacet++)
		{
			pHalfedge = pFacet->halfedge();
			v1 = pHalfedge->vertex()->point();
			v2 = pHalfedge->next()->vertex()->point();
			v3 = pHalfedge->next()->next()->vertex()->point();
			vec.set(v1.x()-sp.x, v1.y()-sp.y, v1.z()-sp.z);
			if (dx.dot(vec)<=0 || dy.dot(vec)<=0 || dz.dot(vec)<=0) {
				continue;
			}
			if (width < fabs(dx.dot(vec))) {
				continue;
			}
			if (height < fabs(dy.dot(vec))) {
				continue;
			}
			vec.set(v2.x()-sp.x, v2.y()-sp.y, v2.z()-sp.z);
			if (dx.dot(vec)<=0 || dy.dot(vec)<=0 || dz.dot(vec)<=0) {
				continue;
			}
			if (width < fabs(dx.dot(vec))) {
				continue;
			}
			if (height < fabs(dy.dot(vec))) {
				continue;
			}
			vec.set(v3.x()-sp.x, v3.y()-sp.y, v3.z()-sp.z);
			if (dx.dot(vec)<=0 || dy.dot(vec)<=0 || dz.dot(vec)<=0) {
				continue;
			}
			if (width < fabs(dx.dot(vec))) {
				continue;
			}
			if (height < fabs(dy.dot(vec))) {
				continue;
			}
			if(bFrontOnly){
				norml.set(pFacet->normal()[0], pFacet->normal()[1], pFacet->normal()[2]);
				if (norml.dot(vec)>viewFrontLow)	continue;
			}
			else if(bSilhouetteOnly){
				norml.set(pFacet->normal()[0], pFacet->normal()[1], pFacet->normal()[2]);
				if (norml.dot(vec)<=viewOscLow || norml.dot(vec)>=viewOscHigh)	continue;
			}
			last_found_Fidx = pFacet->id(); 	
			if(pFacet->have_tag(PICKED)){ //already picked
				continue;
			}
			pFacet->add_tag(PICKED);
			picked_facets.push_back(pFacet->id());
		}
		return last_found_Fidx;
	}
	void unpick_facets_by_column(Vector3 &sp, Vector3 &dx, Vector3 &dy, Vector3 &dz, bool bFrontOnly, bool bSilhouetteOnly=false, FT viewOscLow=-0.0, FT viewOscHigh=1000.0, FT viewFrontLow=0.0)
	{
		float t1 = 0.0f, t2 = 0.0f, t3 = 0.0f;
		float width = dx.magnitude();
		float height = dy.magnitude();
		dx.normalize();
		dy.normalize();
		dz.normalize();
		Facet_iterator pFacet;
		Halfedge_handle pHalfedge;
		Point_3 v1, v2, v3;
		Vector3 vec, vert, norml;
		std::vector<std::list<int>::iterator> it_need_erase;
		for (std::list<int>::iterator it=picked_facets.begin(); it!=picked_facets.end(); it++)
		{
			pFacet = index_to_facet_map[(*it)]; 		
			pHalfedge = pFacet->halfedge();
			v1 = pHalfedge->vertex()->point();
			v2 = pHalfedge->next()->vertex()->point();
			v3 = pHalfedge->next()->next()->vertex()->point();
			vec.set(v1.x()-sp.x, v1.y()-sp.y, v1.z()-sp.z);
			if (dx.dot(vec)<=0 || dy.dot(vec)<=0 || dz.dot(vec)<=0) {
				continue;
			}
			if (width < fabs(dx.dot(vec))) {
				continue;
			}
			if (height < fabs(dy.dot(vec))) {
				continue;
			}
			vec.set(v2.x()-sp.x, v2.y()-sp.y, v2.z()-sp.z);
			if (dx.dot(vec)<=0 || dy.dot(vec)<=0 || dz.dot(vec)<=0) {
				continue;
			}
			if (width < fabs(dx.dot(vec))) {
				continue;
			}
			if (height < fabs(dy.dot(vec))) {
				continue;
			}
			vec.set(v3.x()-sp.x, v3.y()-sp.y, v3.z()-sp.z);
			if (dx.dot(vec)<=0 || dy.dot(vec)<=0 || dz.dot(vec)<=0) {
				continue;
			}
			if (width < fabs(dx.dot(vec))) {
				continue;
			}
			if (height < fabs(dy.dot(vec))) {
				continue;
			}
			if(bFrontOnly){
				norml.set(pFacet->normal()[0], pFacet->normal()[1], pFacet->normal()[2]);
				if (norml.dot(vec)>viewFrontLow)	continue;
			}
			else if(bSilhouetteOnly){
				norml.set(pFacet->normal()[0], pFacet->normal()[1], pFacet->normal()[2]);
				if (norml.dot(vec)<=viewOscLow || norml.dot(vec)>=viewOscHigh)	continue;
			}
			pFacet->remove_tag(PICKED);
			it_need_erase.push_back(it);	
		}
		for (unsigned int i=0; i<it_need_erase.size(); i++)
		{
			picked_facets.erase(it_need_erase[i]);
		}
	}

	void clear_picked_facets(void)
	{
		Facet_handle fh;
		for (std::list<int>::iterator it =picked_facets.begin(); it!=picked_facets.end(); ++it){
			fh = index_to_facet_map[(*it)];
			fh->remove_tag(PICKED);
		}
		picked_facets.clear();
	}

	int nb_picked_facets(void)
	{
		return picked_facets.size();
	}
	void gl_draw_picked_facets(bool smooth_shading, bool use_normals)
	{
		if (picked_facets.size() == 0) {
			return;
		}
		GLfloat mat[4] = {0.85f, 0.6f, 0.1f};
		::glPushAttrib(GL_ALL_ATTRIB_BITS);
		::glDisable(GL_BLEND);
		::glDisable(GL_TEXTURE_2D);
		::glEnable(GL_LIGHTING);
		::glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat);
		::glPolygonMode(GL_FRONT, GL_FILL);		
		::glPolygonMode(GL_BACK, GL_LINE);		
		std::list<int>::iterator pf;
		::glBegin(GL_TRIANGLES);
		for (pf=picked_facets.begin(); pf!=picked_facets.end(); pf++)
			gl_draw_facet(index_to_facet_map[(*pf)], smooth_shading, use_normals, false);
		::glEnd();
		::glPopAttrib();
	}


	FT compute_avg_elen(const std::vector<int> &v)
	{	// compute average edge length in vertex set v
		std::map<int,int> vset;
		std::map<int,int>::iterator iV1, iV2;
		for (int i=0; i<v.size(); i++) {
			vset.insert(std::make_pair(v[i], 0));	// id - bIsHandled
		}
		Halfedge_around_vertex_circulator pHalfEdge, end;
		Vector edge;
		FT len = 0.0;
		int num = 0;
		for (iV1=vset.begin(); iV1!=vset.end(); iV1++)
		{
			pHalfEdge = end = VB(iV1->first);
			CGAL_For_all(pHalfEdge, end)
			{
				if ((iV2=vset.find(pHalfEdge->opposite()->vertex()->id())) != vset.end()) {
					if (iV2->second == 0) {
						const Point& p1 = pHalfEdge->vertex()->point();
						const Point& p2 = pHalfEdge->opposite()->vertex()->point();
						edge = p1 - p2;
						len += (FT)std::sqrt(edge*edge);
						num++;
					}
				}
			}
			iV1->second = 1;
		}
		return (len/(FT)num);
	}

	bool is_boundary_vert(const std::set<int> &pset, int v)
	{
		Halfedge_around_vertex_circulator pHalfEdge = VB(v);
		Halfedge_around_vertex_circulator end = pHalfEdge;
		CGAL_For_all(pHalfEdge, end)
		{
			if (pset.find(pHalfEdge->opposite()->vertex()->id()) == pset.end()) {
				return true;
			}
		}
		return false;
	}

	void sample_vert(const std::vector<int> &vIn, std::vector<int> &vOut, FT dR)
	{
		KD_tree STree;	// kd tree for all vertices
		std::map<Point,int> mVPToIndex;		// vertex point(pos) to index
		std::map<Point,int>::iterator iVPTI;
		for (int i=0; i<vIn.size(); i++) {
			STree.insert(VP(vIn[i]));
			mVPToIndex.insert(std::make_pair(VP(vIn[i]), vIn[i]));
		}
		bool bChange = true;
		std::map<int, short> RM; // indices : removed(1)/processed(2) 
		std::map<int, short>::iterator iRM;
		std::vector<Point> FP;	// found points
		while (bChange)
		{
			bChange = false;
			for (int i=0; i<vIn.size(); i++)
			{
				if (RM.find(vIn[i]) != RM.end())
				{ // found in RM, has been removed or processed
					continue;
				}
				FP.clear();
				Fuzzy_sphere fs(VP(vIn[i]), dR, 0.0);
				STree.search(std::back_inserter(FP), fs);
				for (int j=0; j<FP.size(); j++)
				{	// remove all vertices in fuzzy sphere fs
					if ((iVPTI=mVPToIndex.find(FP[j])) != mVPToIndex.end())
					{	// find the record of FP[i] to get its id
						if (iVPTI->second!=(vIn[i]) && RM.find(iVPTI->second)==RM.end())
						{	// not yet been removed
							RM.insert(std::make_pair(iVPTI->second, 1));	// 1 for remove
							bChange = true;				// changing
						}
					}
				}
				RM.insert(std::make_pair(vIn[i], 2));	// 2 for process
			}
		}
		vOut.clear();
		// collect all the vertices left
		for (int i=0; i<vIn.size(); i++)
		{
			if ((iRM=RM.find(vIn[i])) != RM.end()) {
				if (iRM->second == 1) {	// removed
					continue;
				}
			}
			vOut.push_back(vIn[i]);
		}
	}

	void sample_vert(std::vector<int> &vOut, FT dR)
	{	// from all the vertices of the current mesh
		KD_tree STree;	// kd tree for all vertices
		std::map<Point,int> mVPToIndex;		// vertex point(pos) to index
		std::map<Point,int>::iterator iVPTI;
		for (int i=0; i<size_of_vertices(); i++) {
			STree.insert(VP(i));
			mVPToIndex.insert(std::make_pair(VP(i), i));
		}
		bool bChange = true;
		std::map<int, short> RM; // indices : removed(1)/processed(2) 
		std::map<int, short>::iterator iRM;
		std::vector<Point> FP;	// found points
		while (bChange)
		{
			bChange = false;
			for (int i=0; i<size_of_vertices(); i++)
			{
				if (RM.find(i) != RM.end())
				{ // found in RM, has been removed or processed
					continue;
				}
				FP.clear();
				Fuzzy_sphere fs(VP(i), dR, 0.0);
				STree.search(std::back_inserter(FP), fs);
				for (int j=0; j<FP.size(); j++)
				{	// remove all vertices in fuzzy sphere fs
					if ((iVPTI=mVPToIndex.find(FP[j])) != mVPToIndex.end())
					{	// find the record of FP[i] to get its id
						if (iVPTI->second!=i && RM.find(iVPTI->second)==RM.end())
						{	// not yet been removed
							RM.insert(std::make_pair(iVPTI->second, 1));	// 1 for remove
							bChange = true;				// changing
						}
					}
				}
				RM.insert(std::make_pair(i, 2));	// 2 for process
			}
		}
		vOut.clear();
		// collect all the vertices left
		for (int i=0; i<size_of_vertices(); i++)
		{
			if ((iRM=RM.find(i)) != RM.end()) {
				if (iRM->second == 1) {	// removed
					continue;
				}
			}
			vOut.push_back(i);
		}
	}


	void gl_draw(void)  //only call list
	{
		#if NO_ADV_GPU
			return;
		#endif
		::glPushAttrib(GL_ENABLE_BIT|GL_CURRENT_BIT|GL_HINT_BIT);
		::glCallList(m_gl_display_list);
		::glPopAttrib();
	}

	void gl_draw_face_line(bool smooth_shading, bool use_texcoords)
	{
		::glPushAttrib(GL_ENABLE_BIT|GL_CURRENT_BIT|GL_LINE_BIT|GL_LIGHTING_BIT|GL_POLYGON_BIT);
		// filled polygon
		::glEnable(GL_TEXTURE);
		::glEnable(GL_LIGHTING);
		gl_redraw(smooth_shading, true, use_texcoords);
		// wireframe
		::glDisable(GL_LIGHTING);
		::glDisable(GL_TEXTURE);
		::glColor3f(0.45, 0.45, 0.45);
		::glLineWidth(1.3f);
		::glBegin(GL_LINES);
		for(Facet_iterator pFacet=facets_begin() ; pFacet!=facets_end(); pFacet++)
		{
			Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
			do{
				::glVertex3f(pHalfedge->vertex()->point().x(),pHalfedge->vertex()->point().y(),pHalfedge->vertex()->point().z());
				::glVertex3f(pHalfedge->opposite()->vertex()->point().x(),pHalfedge->opposite()->vertex()->point().y(),pHalfedge->opposite()->vertex()->point().z());
			}
			while( ++pHalfedge != pFacet->facet_begin() );
		}
		::glEnd();
		::glPopAttrib();
	}

	void gl_build_display_list(bool smooth_shading, bool use_normals, bool use_texcoords)
	{
		if (::glIsList(m_gl_display_list)) {
			::glDeleteLists(m_gl_display_list, 1);
		}
		m_gl_display_list = ::glGenLists(1);
		::glNewList(m_gl_display_list, GL_COMPILE);
		gl_redraw(smooth_shading, use_normals, use_texcoords);
		::glEndList();
	}

	// draw	using	OpenGL commands	(display lists)
	void gl_redraw(bool smooth_shading, bool use_normals, bool use_texcoords)
	{
		// draw	polygons
		/*	for(Facet_iterator pFacet=facets_begin(); pFacet!=facets_end(); pFacet++)
		{
		// begin polygon assembly
		::glBegin(GL_POLYGON);
		gl_draw_facet(pFacet,smooth_shading,use_normals);
		::glEnd(); // end polygon assembly
		}*/

		::glBegin(GL_TRIANGLES);
		for(Facet_iterator pFacet=facets_begin(); pFacet!=facets_end(); pFacet++)
		{
			gl_draw_facet(pFacet, smooth_shading, use_normals, use_texcoords);
		}
		::glEnd();
	}

	void gl_draw_facet(Facet_handle pFacet, bool smooth_shading, bool use_normals, bool use_texcoords)
	{
		// one normal per face
		if (use_normals && !smooth_shading)
		{
			const Facet::Normal_3& normal = pFacet->normal();
			glNormal3F(normal[0], normal[1], normal[2]);
		}

		// revolve around current face to get vertices
		Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
		do
		{
			// one normal per vertex
			if (use_normals && smooth_shading)
			{
				const Facet::Normal_3& normal = pHalfedge->vnormal();
				glNormal3F(normal[0], normal[1], normal[2]);
			}
			if (use_texcoords)
			{
				glTexCoord2Fv(pHalfedge->vertex()->texcoord());
			}
			// polygon assembly	is performed per vertex
			const Point& point = pHalfedge->vertex()->point();
			glVertex3F(point[0], point[1], point[2]);
		}
		while(++pHalfedge != pFacet->facet_begin());
	}

	void gl_draw_facet_paramterization (Facet_handle pFacet)
	{
		// revolve around current face to get vertices
		Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
		do
		{
			// polygon assembly	is performed per vertex
			glVertex3F(pHalfedge->vertex()->u(), pHalfedge->vertex()->v(), 0.0);
		}
		while(++pHalfedge != pFacet->facet_begin());
	}

	//////////////////////////////////Begin of VBO///////////////////////////////////
	std::vector<FT> m_pos_array;
	std::vector<FT> m_norm_array;
	void gl_build_vbo(bool smooth_shading, bool use_normals, bool use_texcoords)
	{
		#if NO_ADV_GPU
			return;
		#endif
		if (m_ibo_array == NULL) {
			m_ibo_array = new GLuint [3*size_of_facets()];
		}
		int fi, vi;
		Facet_iterator pFacet;
		tag_vertices(-1);
		m_vb_id.clear();
		m_vb_ia_id.clear();
		for(pFacet=facets_begin(), fi=0; pFacet!=facets_end(); pFacet++, fi++)
		{
			Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
			vi = 0;
			do {
				if (pHalfedge->vertex()->tag()==-1 || pHalfedge->vertex()->dup()==1) {
					const Point& point = pHalfedge->vertex()->point();
					m_pos_array.push_back(point[0]);
					m_pos_array.push_back(point[1]);
					m_pos_array.push_back(point[2]);
					const Facet::Normal_3& normal = pHalfedge->vnormal();
					m_norm_array.push_back(normal[0]);
					m_norm_array.push_back(normal[1]);
					m_norm_array.push_back(normal[2]);
					pHalfedge->vertex()->tag(m_vb_id.size());
					m_vb_id.push_back(pHalfedge->vertex()->id());
					m_vb_ia_id.push_back(3*fi+vi);
				}
				m_ibo_array[3*fi+vi] = pHalfedge->vertex()->tag();
				vi++;
			} while(++pHalfedge != pFacet->facet_begin());
		}

		// Load positions
		glGenBuffers(1, &m_gl_vbo_pos);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		int object_size = sizeof(FT) * 3 * m_vb_id.size();
		glBufferData(GL_ARRAY_BUFFER, object_size, NULL, GL_STATIC_DRAW);
		glBufferSubData(GL_ARRAY_BUFFER, 0, object_size, (FT*)(&(*m_pos_array.begin())));
		// Load normals
		glGenBuffers(1, &m_gl_vbo_norm);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		glBufferData(GL_ARRAY_BUFFER, object_size, NULL, GL_STATIC_DRAW);
		glBufferSubData(GL_ARRAY_BUFFER, 0, object_size, (FT*)(&(*m_norm_array.begin())));
		// Load indices
		glGenBuffers(1, &m_gl_ibo);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo);
		object_size = sizeof(GLuint) * 3 * size_of_facets();
		glBufferData(GL_ELEMENT_ARRAY_BUFFER_ARB, object_size, 0, GL_STATIC_DRAW);
		glBufferSubData(GL_ELEMENT_ARRAY_BUFFER_ARB, 0, object_size, m_ibo_array);
		// End of init
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
		glBindBuffer(GL_ARRAY_BUFFER, 0);
		tag_vertices(0);
	}

	void update_positions_vbo()
	{
#if !NO_ADV_GPU
		if (m_vb_id.empty() || m_vb_ia_id.empty()) {
			return;
		}
		int n_vb_id = m_vb_id.size();
		Point point;
		Facet::Normal_3 normal; 
		for (int i=0; i<n_vb_id; ++i){
			point = VP(m_vb_id[i]);

			m_pos_array[3*i+0] = point[0];
			m_pos_array[3*i+1] = point[1];
			m_pos_array[3*i+2] = point[2];

			normal = VN(m_vb_id[i]);
			m_norm_array[3*i+0] = normal[0];
			m_norm_array[3*i+1] = normal[1];
			m_norm_array[3*i+2] = normal[2];
		}
		//redraw the position buffer
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		int object_size = sizeof(FT) * 3 * n_vb_id;
		glBufferData(GL_ARRAY_BUFFER, object_size, NULL, GL_STATIC_DRAW);
		glBufferSubData(GL_ARRAY_BUFFER, 0, object_size, (FT*)(&(*m_pos_array.begin())));
		//redraw the normal array
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		glBufferData(GL_ARRAY_BUFFER, object_size, NULL, GL_STATIC_DRAW);
		glBufferSubData(GL_ARRAY_BUFFER, 0, object_size, (FT*)(&(*m_norm_array.begin())));
#endif
	}

	void gl_build_vbo_hv(FT *hf)
	{
		if (m_vb_id.empty()) {
			return;
		}
		FT *hv_array = new FT [m_vb_id.size()];
		for (int i=0; i<m_vb_id.size(); i++) {
			hv_array[i] = hf[m_vb_id[i]];
		}
		if (m_gl_vbo_hv == 0) {
			glGenBuffers(1, &m_gl_vbo_hv);
		}
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_hv);
		int object_size = sizeof(FT) * m_vb_id.size();
		glBufferData(GL_ARRAY_BUFFER, object_size, hv_array, GL_STATIC_DRAW);
		// End of init
		glBindBuffer(GL_ARRAY_BUFFER, 0);
		SAFE_DELETE_ARRAY(hv_array);
	}

	void gl_build_vbo_tex_color()
	{
		if (m_vb_id.empty()) {
			return;
		}
		// Build a color buffer for the current VBO
		GLubyte *color_array = new GLubyte [3*m_vb_id.size()];
		for (int i=0; i<m_vb_id.size(); i++) {
			color_array[3*i+0] = m_texture_mat[3*m_vb_id[i]+0];
			color_array[3*i+1] = m_texture_mat[3*m_vb_id[i]+1];
			color_array[3*i+2] = m_texture_mat[3*m_vb_id[i]+2];
		}
		// Load colors
		if (m_gl_vbo_tex_color == 0) {
			glGenBuffers(1, &m_gl_vbo_tex_color);
		}
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_tex_color);
		int object_size = sizeof(GLubyte) * 3 * m_vb_id.size();
		glBufferData(GL_ARRAY_BUFFER, object_size, color_array, GL_STATIC_DRAW);
		// End of init
		glBindBuffer(GL_ARRAY_BUFFER, 0);
		SAFE_DELETE_ARRAY(color_array);
		m_is_built_tex_vbo = true;
	}

	void gl_build_vbo_hf_color(CPseudoColorRGB &PC, FT *hf)
	{
		if (m_vb_id.empty()) {
			return;
		}
		// Build a color buffer for the current VBO
		PC.SetPCValueRange(0.0, 1.0);
		GLubyte *color_array = new GLubyte [3*m_vb_id.size()];
		for (int i=0; i<m_vb_id.size(); i++) {
			PC.GetPC(&color_array[3*i], hf[m_vb_id[i]]);
		}
		// Load colors
		if (m_gl_vbo_hf_color == 0) {
			glGenBuffers(1, &m_gl_vbo_hf_color);
		}
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_hf_color);
		int object_size = sizeof(GLubyte) * 3 * m_vb_id.size();
		glBufferData(GL_ARRAY_BUFFER, object_size, color_array, GL_STATIC_DRAW);
		// End of init
		glBindBuffer(GL_ARRAY_BUFFER, 0);
		SAFE_DELETE_ARRAY(color_array);
	}

	void gl_build_mh_vbo_color(void)
	{
		// Build a color buffer for the current VBO
		GLfloat mat_amb[4];
		GLfloat mat_dif[4];
		GLubyte mat_col[3];
		glGetMaterialfv(GL_FRONT, GL_AMBIENT, mat_amb);
		glGetMaterialfv(GL_FRONT, GL_DIFFUSE, mat_dif);
		mat_col[0] = (mat_amb[0] + mat_dif[0]) * 128;
		mat_col[1] = (mat_amb[1] + mat_dif[1]) * 128;
		mat_col[2] = (mat_amb[2] + mat_dif[2]) * 128;
		GLubyte *color_array = new GLubyte [3*m_vb_id.size()];
		for (int i=0; i<m_vb_id.size(); i++) {
			if (m_hdl[m_vb_id[i]] == -1) {
				color_array[3*i+0] = mat_col[0];
				color_array[3*i+1] = mat_col[1];
				color_array[3*i+2] = mat_col[2];
			} else if (m_hdl[m_vb_id[i]] == m_curr_manip) {
				color_array[3*i+0] = 204;
				color_array[3*i+1] = 51;
				color_array[3*i+2] = 51;
			} else {
				color_array[3*i+0] = 51;
				color_array[3*i+1] = 163;
				color_array[3*i+2] = 51;
			}
		}
		// Load colors
		if (m_gl_vbo_color == 0) {
			glGenBuffers(1, &m_gl_vbo_color);
		}
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_color);
		int object_size = sizeof(GLubyte) * 3 * m_vb_id.size();
		glBufferData(GL_ARRAY_BUFFER, object_size, color_array, GL_STATIC_DRAW);
		// End of init
		glBindBuffer(GL_ARRAY_BUFFER, 0);
		SAFE_DELETE_ARRAY(color_array);
	}

	void gl_draw_hf_vbo(void)
	{
		::glPushAttrib(GL_ENABLE_BIT|GL_LIGHTING_BIT);
		::glEnable(GL_LIGHTING);
		::glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
		::glEnable(GL_COLOR_MATERIAL);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		glVertexPointer(3, GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_VERTEX_ARRAY);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		glNormalPointer(GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_NORMAL_ARRAY);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_hf_color);
		glColorPointer(3, GL_UNSIGNED_BYTE, 0, 0);
		glEnableClientState(GL_COLOR_ARRAY);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo);
		glDrawElements(GL_TRIANGLES, 3*size_of_facets(), GL_UNSIGNED_INT, 0);
		glDisableClientState(GL_NORMAL_ARRAY);
		glDisableClientState(GL_VERTEX_ARRAY);
		glDisableClientState(GL_COLOR_ARRAY);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
		glBindBuffer(GL_ARRAY_BUFFER, 0);
		::glPopAttrib();
	}

	void gl_draw_tex_vbo(void)
	{
		::glPushAttrib(GL_ENABLE_BIT|GL_LIGHTING_BIT);
		::glEnable(GL_LIGHTING);
		::glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
		::glEnable(GL_COLOR_MATERIAL);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		glVertexPointer(3, GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_VERTEX_ARRAY);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		glNormalPointer(GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_NORMAL_ARRAY);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_tex_color);
		glColorPointer(3, GL_UNSIGNED_BYTE, 0, 0);
		glEnableClientState(GL_COLOR_ARRAY);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo);
		glDrawElements(GL_TRIANGLES, 3*size_of_facets(), GL_UNSIGNED_INT, 0);
		glDisableClientState(GL_NORMAL_ARRAY);
		glDisableClientState(GL_VERTEX_ARRAY);
		glDisableClientState(GL_COLOR_ARRAY);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
		glBindBuffer(GL_ARRAY_BUFFER, 0);
		::glPopAttrib();
	}

	void gl_draw_mh_hf_vbo(bool hdl_metaphor)
	{
		::glPushAttrib(GL_ENABLE_BIT|GL_LIGHTING_BIT);
		::glEnable(GL_LIGHTING);
		::glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
		::glEnable(GL_COLOR_MATERIAL);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		glVertexPointer(3, GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_VERTEX_ARRAY);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		glNormalPointer(GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_NORMAL_ARRAY);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_hf_color);
		glColorPointer(3, GL_UNSIGNED_BYTE, 0, 0);
		glEnableClientState(GL_COLOR_ARRAY);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo);
		glDrawElements(GL_TRIANGLES, 3*size_of_facets(), GL_UNSIGNED_INT, 0);
		if (hdl_metaphor)
		{
			glDisableClientState(GL_COLOR_ARRAY);
			for (int i=0; i<m_hdl_f.size(); i++)
			{
				if (i != m_curr_manip) {
					glColor3f(0.2f, 0.64f, 0.2f);
				} else {
					glColor3f(0.8f, 0.2f, 0.2f);
				}
				glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo_mf_hdl[i]);
				glDrawElements(GL_TRIANGLES, 3*m_hdl_f[i].size(), GL_UNSIGNED_INT, 0);
			}
		}
		glDisableClientState(GL_NORMAL_ARRAY);
		glDisableClientState(GL_VERTEX_ARRAY);
		glDisableClientState(GL_COLOR_ARRAY);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
		glBindBuffer(GL_ARRAY_BUFFER, 0);
		::glPopAttrib();
	}

	void gl_draw_vbo(void)
	{
		#if NO_ADV_GPU
			return;
		#endif
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		glVertexPointer(3, GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_VERTEX_ARRAY);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		glNormalPointer(GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_NORMAL_ARRAY);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo);
		glDrawElements(GL_TRIANGLES, 3*size_of_facets(), GL_UNSIGNED_INT, 0);
		glDisableClientState(GL_NORMAL_ARRAY);
		glDisableClientState(GL_VERTEX_ARRAY);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
		glBindBuffer(GL_ARRAY_BUFFER, 0);
	}

	void gl_draw_mh_vbo(bool hdl_metaphor)
	{
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		glVertexPointer(3, GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_VERTEX_ARRAY);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		glNormalPointer(GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_NORMAL_ARRAY);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo);
		glDrawElements(GL_TRIANGLES, 3*size_of_facets(), GL_UNSIGNED_INT, 0);
		if (hdl_metaphor && !m_hdl_f.empty())
		{
			glPushAttrib(GL_ENABLE_BIT|GL_LIGHTING_BIT);
			glEnable(GL_LIGHTING);
			glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE);
			glEnable(GL_COLOR_MATERIAL);
			for (int i=0; i<m_hdl_f.size(); i++)
			{
				if (i != m_curr_manip) {
					glColor3f(0.2f, 0.64f, 0.2f);
				} else {
					glColor3f(0.8f, 0.2f, 0.2f);
				}
				glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo_mf_hdl[i]);
				glDrawElements(GL_TRIANGLES, 3*m_hdl_f[i].size(), GL_UNSIGNED_INT, 0);
			}
			glPopAttrib();
		}
		glDisableClientState(GL_NORMAL_ARRAY);
		glDisableClientState(GL_VERTEX_ARRAY);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
		glBindBuffer(GL_ARRAY_BUFFER, 0);
	}

	void gl_build_hb_vbo_hf_color(CPseudoColorRGB &PC)
	{
		gl_build_vbo_hf_color(PC, m_hb_hf);
	}
	void gl_build_mh_vbo_hf_color(CPseudoColorRGB &PC)
	{
		if (m_curr_manip<0 || m_curr_manip>=m_mh_hf.size()) {
			return;
		}
		gl_build_vbo_hf_color(PC, m_mh_hf[m_curr_manip]);
	}
	void gl_build_hf_s_color(CPseudoColorRGB &PC)
	{
		if (m_is_utd_sb_hf && m_hf_sites.size()>1) {
			gl_build_vbo_hf_color(PC, m_s_hf);
			m_is_built_hf_s_vbo = true;
		}
	}
	void gl_build_hb_vbo_hv(void)
	{
		gl_build_vbo_hv(m_hb_hf);
	}
	void gl_build_mh_vbo_hv(void)
	{
		if (m_curr_manip<0 || m_curr_manip>=m_mh_hf.size()) {
			return;
		}
		gl_build_vbo_hv(m_mh_hf[m_curr_manip]);
	}
	GLuint get_vbo_hv(void)
	{
		return m_gl_vbo_hv;
	}

	void gl_mh_ibo_add(const std::list<Facet_handle> &flist)
	{	// must be called right after 'm_hdl_f.push_back(flist)'
		if (m_ibo_array == NULL) {
			return;
		}
		int num_facets = flist.size();
		GLuint *facet_array = new GLuint [3*num_facets];
		std::list<Facet_handle>::const_iterator pFacet;
		Halfedge_around_facet_circulator pHalfedge;
		int fi;
		for (pFacet=flist.begin(), fi=0; pFacet!=flist.end(); pFacet++, fi++)
		{
			// Revolve around current face to get vertices
			int i = (*pFacet)->facet_begin()->facet()->id();
			facet_array[3*fi+0] = m_ibo_array[3*i+0];
			facet_array[3*fi+1] = m_ibo_array[3*i+1];
			facet_array[3*fi+2] = m_ibo_array[3*i+2];
		}
		// Load indices
		GLuint bufobj;
		glGenBuffers(1, &bufobj);
		m_gl_ibo_mf_hdl.push_back(bufobj);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, bufobj);
		int object_size = sizeof(GLuint) * 3 * num_facets;
		glBufferData(GL_ELEMENT_ARRAY_BUFFER_ARB, object_size, 0, GL_STATIC_DRAW);
		glBufferSubData(GL_ELEMENT_ARRAY_BUFFER_ARB, 0, object_size, facet_array);
		// End of init
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
		// delete data
		SAFE_DELETE_ARRAY(facet_array);
	}
	void gl_mh_ibo_erase(int index)
	{	// must be called right after 'm_hdl_f.erase(m_hdl_f.begin()+m_hdl[(*pv)])'
		glDeleteBuffers(1, &m_gl_ibo_mf_hdl[index]);
		m_gl_ibo_mf_hdl.erase(m_gl_ibo_mf_hdl.begin()+index);
	}

	void gl_draw_vbo_face_line(void)
	{
		#if NO_ADV_GPU
			return;
		#endif
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		glVertexPointer(3, GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_VERTEX_ARRAY);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		glNormalPointer(GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_NORMAL_ARRAY);
		::glPushAttrib(GL_ENABLE_BIT|GL_CURRENT_BIT|GL_LINE_BIT|GL_LIGHTING_BIT|GL_POLYGON_BIT);
		// wireframe
		::glDisable(GL_LIGHTING);
		::glDisable(GL_TEXTURE);
		::glColor3ub(78, 78, 78);
		::glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
		::glLineWidth(1.3f);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo);
		glDrawElements(GL_TRIANGLES, 3*size_of_facets(), GL_UNSIGNED_INT, 0);
		// filled polygon
		::glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
		::glEnable(GL_LIGHTING);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo);
		glDrawElements(GL_TRIANGLES, 3*size_of_facets(), GL_UNSIGNED_INT, 0);
		::glPopAttrib();
	}

	GLfloat m_gl_vbo_line_width;
	GLubyte m_gl_vbo_line_color[3];
	void vbo_line_draw_setting(GLfloat linewidht, GLubyte red,GLubyte green,GLubyte blue){
		m_gl_vbo_line_width = linewidht;
		m_gl_vbo_line_color[0] = red;
		m_gl_vbo_line_color[1] = green;
		m_gl_vbo_line_color[2] = blue;					
	}

	void gl_draw_vbo_line(void)
	{		
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		glVertexPointer(3, GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_VERTEX_ARRAY);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		glNormalPointer(GL_TYPE<FT>::FLAG, 0, 0);
		glEnableClientState(GL_NORMAL_ARRAY);
		::glPushAttrib(GL_ENABLE_BIT|GL_CURRENT_BIT|GL_LINE_BIT|GL_LIGHTING_BIT|GL_POLYGON_BIT);
		// wireframe
		::glDisable(GL_LIGHTING);
		::glDisable(GL_TEXTURE);
		::glColor3ub(m_gl_vbo_line_color[0], m_gl_vbo_line_color[1], m_gl_vbo_line_color[2]);
		::glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
		::glLineWidth(m_gl_vbo_line_width);
		glBindBuffer(GL_ELEMENT_ARRAY_BUFFER_ARB, m_gl_ibo);
		glDrawElements(GL_TRIANGLES, 3*size_of_facets(), GL_UNSIGNED_INT, 0);
		::glPopAttrib();
	}
	void update_vertices_from_bo(GLuint vbo_p, GLuint vbo_n)
	{
		if (m_vb_id.empty() || m_vb_ia_id.empty()) {
			return;
		}
		glBindBuffer(GL_ARRAY_BUFFER, vbo_p);
		float *vp = (float*)glMapBuffer(GL_ARRAY_BUFFER, GL_READ_ONLY);
		if (vp == NULL) {
			glUnmapBuffer(GL_ARRAY_BUFFER);
			return;
		}
		glBindBuffer(GL_ARRAY_BUFFER, vbo_n);
		float *vn = (float*)glMapBuffer(GL_ARRAY_BUFFER, GL_READ_ONLY);
		if (vn == NULL) {
			glUnmapBuffer(GL_ARRAY_BUFFER);
			return;
		}
		FT *pos_array = new FT [3*m_vb_id.size()];
		FT *normal_array = new FT [3*m_vb_id.size()];
		int i0, i1, i2;
		for (int i=0; i<m_vb_id.size(); i++)
		{
			i0 = 3*m_vb_ia_id[i]+0;
			i1 = 3*m_vb_ia_id[i]+1;
			i2 = 3*m_vb_ia_id[i]+2;
			VP(m_vb_id[i]) = Point(vp[i0], vp[i1], vp[i2]);
			pos_array[3*i+0] = vp[i0];
			pos_array[3*i+1] = vp[i1];
			pos_array[3*i+2] = vp[i2];
			normal_array[3*i+0] = vn[i0];
			normal_array[3*i+1] = vn[i1];
			normal_array[3*i+2] = vn[i2];
		}
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		int object_size = sizeof(FT) * 3 * m_vb_id.size();
		glBufferSubData(GL_ARRAY_BUFFER, 0, object_size, pos_array);
		glBindBuffer(GL_ARRAY_BUFFER, vbo_p);
		glUnmapBuffer(GL_ARRAY_BUFFER);
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		object_size = sizeof(FT) * 3 * m_vb_id.size();
		glBufferSubData(GL_ARRAY_BUFFER, 0, object_size, normal_array);
		glBindBuffer(GL_ARRAY_BUFFER, vbo_n);
		glUnmapBuffer(GL_ARRAY_BUFFER);
		glBindBuffer(GL_ARRAY_BUFFER, 0);
		SAFE_DELETE_ARRAY(pos_array);
		SAFE_DELETE_ARRAY(normal_array);
	}

	void update_positions_vbo(GLuint buf_obj)
	{
#if !NO_ADV_GPU
		if (m_vb_id.empty() || m_vb_ia_id.empty()) {
			return;
		}
		glBindBuffer(GL_ARRAY_BUFFER, buf_obj);
		float *vd = (float*)glMapBuffer(GL_ARRAY_BUFFER, GL_READ_ONLY);
		if (vd == NULL) {
			glUnmapBuffer(GL_ARRAY_BUFFER);
			return;
		}
		FT *pos_array = new FT [3*m_vb_id.size()];
		for (int i=0; i<3*size_of_facets(); i++)
		{
			VP(m_vb_id[m_ibo_array[i]]) = Point(vd[3*i+0], vd[3*i+1], vd[3*i+2]);
			pos_array[3*m_ibo_array[i]+0] = vd[3*i+0];
			pos_array[3*m_ibo_array[i]+1] = vd[3*i+1];
			pos_array[3*m_ibo_array[i]+2] = vd[3*i+2];
		}
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_pos);
		int object_size = sizeof(FT) * 3 * m_vb_id.size();
		glBufferSubData(GL_ARRAY_BUFFER, 0, object_size, pos_array);
		glBindBuffer(GL_ARRAY_BUFFER, buf_obj);
		glUnmapBuffer(GL_ARRAY_BUFFER);
		glBindBuffer(0);
		SAFE_DELETE_ARRAY(pos_array);
#endif
	}
	void update_normals_vbo(GLuint buf_obj)
	{
		if (m_vb_id.empty() || m_vb_ia_id.empty()) {
			return;
		}
		glBindBuffer(GL_ARRAY_BUFFER, buf_obj);
		float *vd = (float*)glMapBuffer(GL_ARRAY_BUFFER, GL_READ_ONLY);
		if (vd == NULL) {
			glUnmapBuffer(GL_ARRAY_BUFFER);
			return;
		}
		FT *normal_array = new FT [3*m_vb_id.size()];
		for (int i=0; i<3*size_of_facets(); i++)
		{
			normal_array[3*m_ibo_array[i]+0] = vd[3*i+0];
			normal_array[3*m_ibo_array[i]+1] = vd[3*i+1];
			normal_array[3*m_ibo_array[i]+2] = vd[3*i+2];
		}
		glBindBuffer(GL_ARRAY_BUFFER, m_gl_vbo_norm);
		int object_size = sizeof(FT) * 3 * m_vb_id.size();
		glBufferSubData(GL_ARRAY_BUFFER, 0, object_size, normal_array);
		glBindBuffer(GL_ARRAY_BUFFER, buf_obj);
		glUnmapBuffer(GL_ARRAY_BUFFER);
		glBindBuffer(0);
		SAFE_DELETE_ARRAY(normal_array);
	}
	//////////////////////////////////End of VBO/////////////////////////////////////

	// superimpose edges
	void superimpose_edges(bool	skip_ordinary_edges=true,
		bool skip_control_edges = false,
		bool voronoi_edge = false)
	{
		::glBegin(GL_LINES);
		for (Edge_iterator h=edges_begin(); h!=edges_end(); h++)
		{
			// ignore	this edges
			if(skip_ordinary_edges && !h->control_edge())
				continue;

			// ignore	control	edges
			if(skip_control_edges && h->control_edge())
				continue;

			if(voronoi_edge)
			{
				Facet_handle pFace1 = h->facet();
				Facet_handle pFace2 = h->opposite()->facet();
				if(pFace1 == NULL || pFace2 == NULL)
					continue;

				const Point &p1 = h->vertex()->point();
				const Point &p2 = h->next()->vertex()->point();
				const Point &p3 = h->next()->next()->vertex()->point();

				kernel k;
				Point d1 = k.construct_circumcenter_3_object()(p1,p2,p3);
				glVertex3F(d1[0],d1[1],d1[2]);

				const Point &pp1 = h->opposite()->vertex()->point();
				const Point &pp2 = h->opposite()->next()->vertex()->point();
				const Point &pp3 = h->opposite()->next()->next()->vertex()->point();
				Point d2 = k.construct_circumcenter_3_object()(pp1,pp2,pp3);
				glVertex3F(d2[0],d2[1],d2[2]);
			}
			else
			{
				// assembly	and	draw line	segment
				const Point& p1 = h->prev()->vertex()->point();
				const Point& p2 = h->vertex()->point();
				glVertex3F(p1[0],p1[1],p1[2]);
				glVertex3F(p2[0],p2[1],p2[2]);
			}
		}
		::glEnd();
	}

	// superimpose vertices
	void superimpose_vertices()
	{
		::glBegin(GL_POINTS);
		for(Point_iterator pPoint = points_begin();
			pPoint !=	points_end();
			pPoint++)
			glVertex3F(pPoint->x(),pPoint->y(),pPoint->z());
		::glEnd(); //	// end point assembly
	}

	// superimpose vertices
	void superimpose_spheres(FT scale)
	{
		GLUquadricObj* pQuadric = gluNewQuadric();
		::glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
		for (Vertex_iterator pVertex=vertices_begin(); pVertex!=vertices_end(); pVertex++)
		{
			::glPushMatrix();
			FT radius = average_edge_length_around(pVertex);
			::glTranslated(pVertex->point().x(), pVertex->point().y(), pVertex->point().z());
			::gluSphere(pQuadric,scale*radius,24,24); 
			::glPopMatrix();
		}
		gluDeleteQuadric(pQuadric);
	}

	// write in	obj	file format	(OBJ).
	void write_obj(char	*pFilename, int incr=1) //	1-based	by default
	{
		if (pFilename == NULL) {
			return;
		}
		std::ofstream	stream(pFilename);
		locale NarrowCharStyle("C");
		stream.imbue(NarrowCharStyle);
		// output	vertices
		for (Point_iterator pPoint=points_begin(); pPoint!=points_end(); pPoint++) {
			stream <<'v'<<' '<<pPoint->x()<<' '<<pPoint->y()<<' '<<pPoint->z()<< std::endl;
		}

		// precompute	vertex indices
		if (!m_is_set_items_id) {
			set_hds_items_id();
		}

		// output	facets
		for(Facet_iterator pFacet	=	facets_begin();
			pFacet !=	facets_end();	
			pFacet++)	
		{
			stream <<	"f";
			Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
			do 
			stream <<	" "	<< pHalfedge->vertex()->id()+incr;
			while(++pHalfedge	!= pFacet->facet_begin());
			stream <<	std::endl;
		}
		stream.close();
	}

	// write in	ply2 file format
	void write_ply2(char *pFilename) //	1-based	by default
	{
		if (pFilename == NULL) {
			return;
		}
		std::ofstream stream(pFilename);

		stream <<size_of_vertices()<< std::endl <<size_of_facets()<< std::endl;
		// output vertices
		for (Point_iterator pPoint=points_begin(); pPoint!=points_end(); pPoint++) {
			stream <<pPoint->x()<<' '<<pPoint->y()<<' '<<pPoint->z()<< std::endl;
		}
		// precompute vertex indices
		if (!m_is_set_items_id) {
			set_hds_items_id();
		}
		// output facets
		for (Facet_iterator pFacet=facets_begin(); pFacet!=facets_end(); pFacet++)	
		{
			Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
			stream << CGAL::circulator_size(pHalfedge) << ' ' << pHalfedge->vertex()->id();
			while(++pHalfedge != pFacet->facet_begin()) {
				stream <<' '<<pHalfedge->vertex()->id();
			}
			stream << std::endl;
		}
	}



	// count #boundaries
	unsigned int nb_boundaries()
	{
		unsigned int nb = 0;
		tag_halfedges(0);
		Halfedge_handle	seed_halfedge	=	NULL;
		while((seed_halfedge = get_border_halfedge_tag(0)) !=	NULL)
		{
			nb++;
			seed_halfedge->tag(1);
			Vertex_handle	seed_vertex	= seed_halfedge->prev()->vertex();
			Halfedge_handle	current_halfedge = seed_halfedge;
			Halfedge_handle	next_halfedge;
			do
			{
				next_halfedge	=	current_halfedge->next();
				next_halfedge->tag(1);
				current_halfedge = next_halfedge;
			}
			while(next_halfedge->prev()->vertex()	!= seed_vertex);
		}
		return nb;
	}

	// get any border	halfedge with	tag
	Halfedge_handle	get_border_halfedge_tag(int	tag)
	{
		for(Halfedge_iterator pHalfedge	=	halfedges_begin();
			pHalfedge	!= halfedges_end();
			pHalfedge++)
			if(pHalfedge->is_border()	&&
				pHalfedge->tag()	== tag)
				return pHalfedge;
		return NULL;
	}

	// get any facet with tag
	Facet_handle get_facet_tag(const int tag)
	{
		for (Facet_iterator pFace = facets_begin(); pFace != facets_end(); pFace++)
			if (pFace->tag() == tag)
				return pFace;
		return NULL;
	}

	// tag component 
	void tag_component(Facet_handle	pSeedFacet, const int tag_free, const int tag_done)
	{
		pSeedFacet->tag(tag_done);
		std::list<Facet_handle> facets;
		facets.push_front(pSeedFacet);
		while(!facets.empty())
		{
			Facet_handle pFacet = facets.front();
			facets.pop_front();
			pFacet->tag(tag_done);
			Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
			Halfedge_around_facet_circulator end = pHalfedge;
			CGAL_For_all(pHalfedge,end)
			{
				Facet_handle pNFacet = pHalfedge->opposite()->facet();
				if(pNFacet != NULL && pNFacet->tag() == tag_free)
				{
					facets.push_front(pNFacet);
					pNFacet->tag(tag_done);
				}
			}
		}
	}

	// count #components
	unsigned int nb_components()
	{
		unsigned int nb = 0;
		tag_facets(0);
		Facet_handle seed_facet = NULL;
		while((seed_facet = get_facet_tag(0)) != NULL)
		{
			nb++;
			tag_component(seed_facet,0,1);
		}
		return nb;
	}

	// compute the genus
	// V - E + F + B = 2 (C	- G)
	// C ->	#connected components
	// G : genus
	// B : #boundaries
	int	genus()
	{
		int	c	=	nb_components();
		int	b	=	nb_boundaries();
		int	v	=	size_of_vertices();
		int	e	=	size_of_halfedges()/2;
		int	f	=	size_of_facets();
		return genus(c,v,f,e,b);
	}
	int	genus(int c, int v, int f, int e, int b)
	{
		return (2*c+e-b-f-v)/2;
	}
	std::vector< std::vector<Point> > m_regulateV;
	bool m_bShowRegulateTempV;
	void facetRegulate_init(){
		m_bShowRegulateTempV = false;
		int nf = size_of_facets();
		m_regulateV.resize(nf);
		for (int i=0; i<nf; i++)	m_regulateV[i].resize(3);
	}
	void facetRegulate(){
		int fid = 0;
		for (Facet_iterator pFacet=facets_begin(); pFacet!=facets_end(); pFacet++)	
		{
			Halfedge_around_facet_circulator h = pFacet->facet_begin();
			int vi = 0;
			CGAL_For_all(h, pFacet->facet_begin() ){
				Vector edge = (h->next()->vertex()->point() - h->vertex()->point());
				FT eLen = std::sqrt(edge*edge);
				FT offset = eLen * 0.5 / std::sqrt(3.0f);
				Point edgeCenter = h->vertex()->point()+ edge * 0.5; 
				Vector outENormal = CGAL::cross_product(edge, pFacet->normal());
				FT oenLen = std::sqrt(outENormal*outENormal);
				outENormal = outENormal / oenLen * offset;
				Point np = edgeCenter + outENormal;
				m_regulateV[fid][vi] = np;
				vi++;
			}
			Vector3 newCenter = ( P2V3(m_regulateV[fid][0]) + P2V3(m_regulateV[fid][1]) + P2V3(m_regulateV[fid][2]) ) /3;
			Point vv[3];
			for (int i=0; i<3; i++) vv[i] = m_regulateV[fid][i];
			m_regulateV[fid][0] = V32P( newCenter + (newCenter-P2V3(vv[1])) );
			m_regulateV[fid][1] = V32P( newCenter + (newCenter-P2V3(vv[2])) );
			m_regulateV[fid][2] = V32P( newCenter + (newCenter-P2V3(vv[0])) );
			fid++;
		}
	}
	void stitching_regulate(){
		int vid = 0;
		for (Vertex_iterator pVertex=vertices_begin(); pVertex!=vertices_end(); pVertex++, vid++)	
		{
			Halfedge_around_vertex_circulator he = pVertex->vertex_begin();
			int nbde = 0;
			Vector3 newPoint(0,0,0);
			CGAL_For_all(he, pVertex->vertex_begin()){
				if(!he->is_border()){
					int fid = he->facet()->id();
					for (int i=0; i<3; i++)
					{
						if( fvid(fid, i) == vid )	{
							newPoint += P2V3 (m_regulateV[fid][i] );
							break;
						}
					}	
				}
			}
			Point& point = pVertex->point();
			int nf = is_border(pVertex)? valence(pVertex)-1 : valence(pVertex);
			point = V32P( newPoint/float(nf) );
		}
		compute_normals();
		update_positions_vbo();
	}

	void gl_draw_napolean(CPseudoColorRGB &PC)
	{
		GLubyte PsdColor[3]= {255, 255, 255};
		PC.SetPCValueRange(0.0, 1.0);
		FT abs_max = m_regulateV.size();
		::glPointSize(5.5);
		if (m_regulateV.size() == 0) {	return;	}
		::glPushAttrib(GL_ALL_ATTRIB_BITS/*GL_LIGHTING_BIT|GL_CURRENT_BIT*/);
		::glDisable(GL_LIGHTING);
		::glDisable(GL_TEXTURE_2D);
		::glBegin(GL_POINTS);	
		for (int i=0; i<m_regulateV.size(); ++i) {
			PC.GetPC(PsdColor, (FT)i/(FT)abs_max);
			::glColor3ubv(PsdColor);
			for (int j=0; j<m_regulateV[i].size(); j++)
			{
				glVertex3f(m_regulateV[i][j].x(), m_regulateV[i][j].y(), m_regulateV[i][j].z());
			}
		}
		::glEnd();
		::glPopAttrib();
	}
	void maxError_regulate(){
		compute_all_edge_len();
		FT maxDeviation = 0;
		for(Facet_iterator fit = facets_begin(); fit != facets_end(); ++fit){

			Halfedge_around_facet_circulator he = fit->facet_begin();
			std::vector<FT> el;
			CGAL_For_all(he, fit->facet_begin()){
				el.push_back( m_edge_Len[he->id()/2] );  
			}
			std::deque<FT>::iterator it1, it2;;
			FT maxEdiff = 0;
			FT totalELen = 0;
			for (int j=0; j<3; j++)
			{
				if( fabs( el[j]-el[(j+1)%3] ) > maxEdiff )
					maxEdiff = fabs(el[j] - el[(j+1)%3]);
				totalELen += el[j];
			}
			FT avgElen = totalELen/3.0f;
			FT deviation = maxEdiff / avgElen ;
			if( deviation > maxDeviation )	maxDeviation = deviation;
		}
		std::cerr<<maxDeviation<<"\n";
	}
};

template <class	kernel, class items>
class ROI_vertex
{
	typedef typename Enriched_polyhedron<kernel,items>::Vertex_handle Vertex_handle;
private:
	Vertex_handle	m_v;
	bool			m_border;

public:
	ROI_vertex()
	{
		m_v = 0;
		m_border = false;
	}
	ROI_vertex(Vertex_handle &v, bool b=false)
	{
		m_v = v;
		m_border = b;
	}
	Vertex_handle&	vertex() {return m_v;}
	void			set_vertex_handle(Vertex_handle &v) {m_v = v;}
	bool			is_border() {return m_border;}
	void			set_border(bool b) {m_border = b;}
};

// compute facet normal
template <class	Kernel, class Facet>
struct Facet_normal	// (functor)
{
	typedef typename Kernel::FT FT;
	typedef typename typename Facet::Halfedge_handle Halfedge_handle;
	typedef typename Facet::Normal_3 FNormal;
	void operator()(Facet& f)
	{
		FNormal sum = CGAL::NULL_VECTOR;
		Facet::Halfedge_around_facet_circulator h = f.facet_begin();
		do {
			FNormal normal	= CGAL::cross_product(
				h->next()->vertex()->point() - h->vertex()->point(),
				h->next()->next()->vertex()->point() - h->next()->vertex()->point());
			FT sqnorm = normal * normal;
			if(sqnorm != 0)
				normal = normal	/ (float)std::sqrt(sqnorm);
			sum = sum + normal;
		} while (++h != f.facet_begin());
		FT sqnorm = sum * sum;
		if (sqnorm != 0.0) {
			f.normal() = sum / std::sqrt(sqnorm);
		} else {
			f.normal() = CGAL::NULL_VECTOR;
			TRACE("degenerate face\n");
		}
	}
};

// compute facet normal	
template <class	Kernel, class Facet>
struct Facet_normal2	// (functor)
{
	typedef typename Kernel::FT FT;
	typedef typename typename Facet::Halfedge_handle Halfedge_handle;
	typedef typename Facet::Normal_3 FNormal;
	void operator()(Facet& f)
	{
		Halfedge_handle h = f.halfedge();
		if (h == NULL) {
			f.normal() = CGAL::NULL_VECTOR;
			return;
		}
		FNormal normal	= CGAL::cross_product(
			h->next()->vertex()->point() - h->vertex()->point(),
			h->next()->next()->vertex()->point() - h->next()->vertex()->point());
		FT sqnorm = normal * normal;
		if (sqnorm != 0.0) {
			f.normal() = normal	/ (float)std::sqrt(sqnorm);
		} else {
			f.normal() = CGAL::NULL_VECTOR;
			TRACE("degenerate face\n");
		}
	}
};

// compute vertex normal 
template <class	Kernel, class Vertex>
struct Vertex_normal //	(functor)
{
	typedef typename Kernel::FT FT;
	typedef typename Vertex::Normal Normal;
	void operator()(Vertex&	v)
	{
		Normal	normal = CGAL::NULL_VECTOR;
		Vertex::Halfedge_around_vertex_const_circulator	pHalfedge =	v.vertex_begin();
		Vertex::Halfedge_around_vertex_const_circulator	begin = pHalfedge;
		CGAL_For_all(pHalfedge,begin) {
			if(!pHalfedge->is_border()) {
				normal = normal	+ pHalfedge->facet()->normal();
			}
		}
		FT sqnorm = normal * normal;
		if (sqnorm != 0.0f) {
			v.normal() = normal	/ (float)std::sqrt(sqnorm);
		} else {
			v.normal() = CGAL::NULL_VECTOR;
		}
	}
};

// compute vertex normal for rendering
template <class	Kernel, class Vertex>
class HalfedgeVertex_normal //	(functor)
{
	typedef typename Kernel::FT FT;
	typedef typename Vertex::Halfedge_around_vertex_circulator HCirculator;
	typedef typename Vertex::Normal Normal;
private:
	FT m_threshold;
public:
	HalfedgeVertex_normal(FT angle) : m_threshold(Cos(angle)) {}
	void operator()(Vertex& v)
	{
		HCirculator pHalfedge =	v.vertex_begin();
		if (pHalfedge == NULL) {	// isolated vertex
			return;
		}
		HCirculator begin = pHalfedge;
		std::vector<HCirculator> ridge;
		FT max_cos = m_threshold;
		HCirculator pHEMax;
		Normal normal;
		// find ridge edges
		CGAL_For_all(pHalfedge, begin)
		{
			if (pHalfedge->is_border() || pHalfedge->opposite()->is_border()) {
				ridge.push_back(pHalfedge);
				continue;
			}
			FT cos_value = pHalfedge->facet()->normal() * pHalfedge->opposite()->facet()->normal();
			if (cos_value < m_threshold) {
				ridge.push_back(pHalfedge);
			} else if (cos_value >= max_cos) {
				max_cos = cos_value;
				pHEMax = pHalfedge;
			}
		}
		if (ridge.size() == 1) {
			ridge.push_back(pHEMax);
		}
		// no ridge edge, compute average of all faces
		if (ridge.size() == 0) {
			//normal = CGAL::NULL_VECTOR;
			//pHalfedge =	v.vertex_begin();
			//CGAL_For_all(pHalfedge, begin) {
			//	if(!pHalfedge->is_border()) {
			//		normal = normal + pHalfedge->facet()->normal();
			//	}
			//}
			//FT sqnorm = normal * normal;
			//if (sqnorm != 0.0f) {
			//	normal = normal / (FT)std::sqrt(sqnorm);
			//} else {
			//	normal = begin->facet()->normal();
			//}
			pHalfedge =	v.vertex_begin();
			CGAL_For_all(pHalfedge, begin) {
				if(!pHalfedge->is_border()) {
					//pHalfedge->vnormal() = normal;
					pHalfedge->vnormal() = v.normal();	// average normal has been computed
				}
			}
			v.dup(0);	// no ridge, no need to duplicate
			return;
		}
		// compute average in each group (divided by ridges)
		for (int i=0; i<static_cast<int>(ridge.size()); i++)
		{
			if (ridge[i]->is_border()) {
				continue;
			}
			normal = CGAL::NULL_VECTOR;
			pHalfedge =	ridge[i];
			CGAL_For_all(pHalfedge, ridge[(i+1)%ridge.size()]) {
				normal = normal + pHalfedge->facet()->normal();
			}
			FT sqnorm = normal * normal;
			if (sqnorm != 0.0f) {
				normal = normal / (FT)std::sqrt(sqnorm);
			} else {
				normal = ridge[i]->facet()->normal();
			}
			pHalfedge =	ridge[i];
			CGAL_For_all(pHalfedge, ridge[(i+1)%ridge.size()]) {
				pHalfedge->vnormal() = normal;
			}
		}
		v.dup(1);	// need to duplicate
	}
};

